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Abstract 


The  high  temperature  (1500  n^K)  exhaust  gases  from  an  airborne 
chemical  laser  being  forcibly  ejected  downwards  at  a  jet  to  freestream 
dynamic  pressure  ratio  (Q)  of  0.15  from  an  aspect  ratio  1*75  rectangular 
diffuser  orifice  aligned  with  the  major  axis  parallel  to  the  Mach  0.7 
ambient  crossflow  at  the  tropopause  was  successfully  simulated  using  a 
modified  donor-cell  ICE  (implicit ^Continuous  fluid  jhxlerian)  algorithm 
formulated  in  primitive  variables.  The  complete  set  of  time  dependent, 
three-dimensional  Navier-Stokes  equations  and  a  species  conservation 
'  equation  were  numerically  solved  for  every  cell  in  the  computational 
domain.  The  diffusive  flux  effects  caused  by  concentration  gradients 
(Fick's  Law)  as  well  as  variable  transport  and  thermodynamic  properties 
of  the  gas  mixture  were  incorporated  into  the  numerical  model. 
Turbulence  closure  was  achieved  by  employing  a  locally  varying  velocity 
defect  eddy  viscosity  model.  Chemical  reactions  between  the  exhaust 
gases  and  the  ambient  crossflow  were  proscribed.^ 

''"'Data  acquired  from  the  numerical  simulations  were  used  to  define 
the  trajectory  of  the  jet  plume,  the  extent  of  the  recirculation  zone  in 
the  wake,  and  the  regions  with  possibly  large  heat  transfer  rates 
Simplified  analyses  were  also  conducted  to  determine  whether  essential 
flow  phenomena  were  captured  mathematically  using  simple  binary  gas 
interactions.  Convective  processes  were  observed  to  dominate  the  low  Q 
jet-crossflow  interaction  for  gas  mixture,  helium,  and  air  injectants; 
each  of  the  simulations  contained  essentially  the  same  flow  features  and 
characteristics.  Thermal  diffusion  was  also  seen  to  have  a 


significantly  greater  effect  than  molecular  diffusion  for  the  jet- 


xv 


AFIT/DS/AA/80-1 


crossflow  gases  simulated.  For  the  low  molecular  weight  injectants,  a 
secondary  wall  vortex  (recirculation  2one)  was  created  approximately  3.2 
jet  streamwise  dimensions  (D)  aft  of  the  jet  center  apparently  caused  by 
diffusive  transport  due  to  concentration  gradients  in  the  presence  of  a 
constraining  "no-slip"  wall.  The  lateral  span  of  che  wall  vortex  did 
not  exceed  0.3  D. 

Jet  penetration  of  the  ambient  flowfield  was  observed  to  be 
dependent  upon  the  molecular  weight  of  the  injectant  for  the  constant  Q 
constraint.  A  molecular  weight  correction  factor  used  in  typical  heat- 
mass  transfer  boundary  layer  problems  was  utilized  to  synthesize  the 
trajectory  curve  of  one  gas  from  that  of  another  gas  at  the  same 
conditions  and  to  correct  empirically  derived  trajectory  formulae  for 
variances  in  molecular  weight.  The  latter  technique  resulted  in  less 
than  5  percent  variation  from  the  corresponding  trajectory  derived  from 
the  numerical  simulation.  Sensitivity  analyses  relating  the  heat 
transfer  to  the  injection  surface  from  the  jet  plume  with  the  magnitude 
of  the  turbulent  dif fusivities  were  also  conducted. 


TIME  DEPENDENT  NAVIER-STOKES  SOLUTION 
OF  A  TURBULENT  GAS  JET  EJECTED  FROM  A  RECTANGULAR 
ORIFICE  INTO  A  HIGH-SUBSONIC  CROSS  FLOW 

I.  Introduction 

The  development  of  the  Airborne  leaser  Laboratory  II  (  ALL  II  ) 
involves  the  emplacement  of  a  large  chemical  laser  on  board  a  wide  body 
jet  aircraft.  This  study  in  support  of  the  ALL  II  program  is  concerned 
with  the  three-dimensional  flowfield  produced  by  the  interaction  of  the 
high  temperature,  low  molecular  weight  waste  gas  (generated  by  the 
chemical  laser  system  and  forcibly  exhausted  normally  downward  from  a 
rectangular  diffuser  embedded  in  the  aft  fuselage)  with  the  surrounding 
high-subsonic  ambient  freestream.  Because  of  both  the  high  temperature 
and  the  corrosive  nature  of  the  chemical  laser  exhaust,  it  is  imperative 
to  determine  the  extent  of  the  heating  on  the  aircraft  -uselage  and  the 
plume  of  the  laser  effluent.  Consequently,  the  primary  goals  of  this 
research  effort  are  (1)  to  demonstrate  the  feasibility  of  the  numerical 
methodology  to  accurately  predict  the  viscous  jet-crossflow  interaction, 
(2)  to  solve  the  three-dimensional  jet  trajectory  problem  for  the  real 
gas  exhaust  mixture,  and,  (3)  to  determine  the  location  and  extent  of 
regions  with  possibly  large  rates  of  heat  transfer  to  the  injection 
surface.  Ancillary  goals  of  determining  the  location  and  extent  of  any 
recirculation  zone,  the  structure  and  characteristics  of  the  jet- 
crossflow  interaction,  and  the  effects  of  the  diffusivity  induced  by  the 
turbulence  on  the  heat  transfer  from  the  hot  gas  mixture  co  the  aircraft 
surface  are  also  of  prime  importance. 

It  is  widely  recognized  that  there  exist  at  least  three  principal 
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regions  of  interest  for  problems  in  which  a  turbulent  jet  is  injected 
into  a  crossflow  independent  of  the  respective  gases  as  indicated  in 
Figure  1.  The  first  of  these  is  designated  as  the  "near  field"  or 
"potential  core"  region  and  is  characterized  by  the  irrotationality  and 
dominance  of  the  jet  flow.  As  the  crossflow  impacts  the  issuing  jet,  a 
compression  zone  is  formed  ahead  of  the  jet  while  a  pressure  defect  zone 
is  formed  behind  the  jet.  It  is  this  induced  pressure  difference  which 
initiates  the  deformation  of  the  jet  (Ref  1)  and  which  consequently 
provides  the  basis  for  the  interactive  mixing  of  the  two  flows.  The 
second,  or  intermediate,  region  has  been  termed  "curvilinear"  by  Chan 
and  Kennedy  (Ref  2)  and  "zone  of  maximum  deflection"  by  Fratte  and 
Baines  (Ref  3) •  It  is  characterized  by  the  rapid  change  in  direction  of 
the  jet  from  its  initial  injection  to  that  of  the  governing  crossflow. 
Within  this  zone,  a  pair  of  counterrotating  vortices  are  formed  which 
enhance  the  entrainment  of  fluid  from  the  crossflow  into  the  jet  plume 
and  which  directly  affect  the  jet-crossflow  interaction.  The  final 
region  is  designated  as  "far  field"  in  which  the  jet  flow  is  essentially 
aligned  with  the  crossflow  and  in  which  it  Is  assumed  that  profiles  of 
velocity  and  fluid  properties  exhibit  Gaussian  similitude;  additionally, 
the  flow  is  further  assumed  to  be  isotropic. 

Real  gases  at  elevated  temperatures  being  forcibly  ejected  into  a 
cross  flowing  stream  of  another  gas  complicate  any  mathematical  modeling 
of  the  flow  interaction.  Adequate  description  of  the  gas  dynamics  is 
proscribed  by  approximate  integral,  potential,  or  boundary  layer 
methodologies;  consequently,  only  the  full  set  of  time  dependent  Navier- 
Stokes  equations  utilizing  both  variable  gas  transport  and  thermodynamic 
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properties  become  requisite.  Prior  to  this  study  and  the  advancement  in 
technology  it  represents,  no  documented  numerical  method  existed  which 
could  provide  time  dependent  solutions  to  problems  of  this  type  in 
three-dimensions.  The  development  and  validation  of  this  numerical 
scheme  as  a  useful  engineering  predictive  tool  within  each  of  the 
regions  associated  with  the  jet-crossflow  interaction  comprise  the  bulk 
of  the  research  effort. 

Historically,  the  early  studies  of  jets  injected  into  crossflows 
were  concerned  with  obtaining  empirical  formulae  to  define  the  jet 
trajectory.  There  exist  many  means  of  mathematically  defining  the 
trajectory  curve,  including  those  of  the  following  list: 

(1)  Locus  of  maximum  velocity, 

(2)  Locus  of  maximum  total  pressure, 

C3)  Locus  of  maximum  total  temperature, 

(4)  Locus  of  maximum  species  concentration, 

(5)  Locus  of  maximum  /  minimum  density, 

(6)  Locus  of  maximum  pressure  coefficient. 

Host  researchers  have  normally  chosen  to  define  the  mean  flow  trajectory 
by  curve  fitting  the  loci  of  maximum  velocity  or  momentum  flux. 
Additionally,  the  empirical  methodology  is  further  split  into  two 
groups:  (1)  analytical,  and  (2)  experimental. 

The  major  impetus  affecting  jet-crossflow  research  was  the  advent 
of  modern  V/STOL  aircraft  and  the  attendant  need  to  define  the 
aerodynamic  impact  of  the  interaction  of  the  jet-crossflow  upon  the  air 
vehicle  performance,  stability  and  control.  Not  only  ware  trajectory 
data  of  prime  importance,  but  also  data  describing  the  pressure  and 
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velocity  fields  near  the  jet  orifice  and  their  associated  effects  on  the 
vehicle  surface  became  equally  vital.  Paralleling  the  advance  in 
aircraft  technology,  large  commercial  power  stations  came  into  vogue  and 
the  need  to  define  the  impact  of  exhaust  particulate  dispersion  on  the 
local  environment  also  became  a  matter  of  research  interest  as  well  as 
civic  concern.  These  two  driving  forces  helped  focus  research  effort  on 
approximate  integral  and  potential  flow  methodologies.  Finally,  the 
Navier-Stokes  solutions  of  the  jet-crossflow  problem  became  more 
tractable  with  the  development  of  modern  computational  machinery.  Each 
of  the  aforementioned  methodologies  will  be  briefly  examined  in  the  next 
section. 

1.1  Jets  in  Crossflows  -  Background 

Experiment  and  Empirical  Methods.  The  early  analytical  studies 
detailing  the  behavior  of  jet  injection  into  a  crossflow  were 
accomplished  for  round  jets  because  of  the  reduction  in  complexity 
afforded  by  the  geometry.  The  simplest  analyses  were  concerned  with 
deriving  mathematical  expressions  for  the  jet  trajectory.  Abroraovich 
used  the  method  developed  by  M.  S.  Volinskii  (Ref  1)  in  which  each 
infinitesimal  segment  of  the  jet  is  treated  as  an  airfoil  at  some  given 
sweep  angle  with  respect  to  the  crossflow  (see  Figure  2.).  Enforcing 
the  condition  of  radial  equilibrium,  the  normal  component  of  the 
sectional  aerodynamic  drag  was  equated  to  the  centrifugal  force  from 
which  the  local  radius  of  curvature  of  the  jet  trajectory  was  computed. 
Abromovich's  following  expression  (12.162)  for  the  trajectory  of  a  jet 
normally  injected  from  a  flat  plate  into  a  crossflow  results: 
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d  is  the  initial  diameter  of  the  jet  at 
injection; 

x,y  are  the  streamwi&e  and  normal  coordinate 
directions,  respectively; 


and  is  the  jet  element  drag  coefficient. 

Vizel  and  Mostinskii  (Ref  2)  improved  the  agreement  of  this  semi- 
empirical  method  by  determining  the  trajectory  of  the  jet  utilizing  the 
drag  coefficient  of  the  entire  jet,  C^.  If  experimental  results  were 
available,  the  following  expression 
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could  be  used  to  determine  the  trajectory  of  the  jet.  Further 

refinements  to  this  methodology  were  added  by  Epshtein  (Ref  5)  and 
Vakhlamov  (Ref  6).  Epshtein  removed  the  restriction  of  the  normal 

component  of  the  jet  momentum  being  constant  and  permitted  it  to  vary 
with  the  gravitational  force.  Vakhlamov  developed  prescribed  control 


surfaces  over  which  a  momentum  balance  was  conducted.  His  resulting 
expression  for  the  trajectory  of  a  round  jet  injected  normally  into  a 
crossflow  is  given  by 
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where  S  ■  x/d  and  9  “  y/d. 


Shandarov  (Ref  7)  following  the  basic  methodology  used  by  Abromovlch, 
obtained  a  similar  expression  for  circular  jets  ejected  normally  into  a 
uniform  unbounded  subsonic  gas  flow, 
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where  q^  *  pUoo  and  q2  =“  Pj  Vj  at  injection. 

Despite  the  degree  of  analysis  each  of  the  derivations  for  the 
preceding  expressions  entailed,  the  range  of  jet-crossflow  parameters  is 
rather  limited  to  be  of  practical  engineering  value.  Consequently, 
empirical  power  )aw  expressions  for  the  trajectory  have  been  developed 
for  a  wider  range  of  temperature,  density,  and  velocity  ratios. 
Abromovich  presented  two  empirical  expressions  developed  by  Shandarov 


and  Ivanov  ostensibly  for  air  jets.  They  are  given  by  the  following 
expressions  for  normal  jet-crosaf low  injection: 


Shandarov  : 
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Callaghan  and  Ruggeri  (Ref  8)  experimentally  determined  an  empirical 
expression  valid  for  heated  (860  °R)  circular  air  jets  with  jet  Reynolds 
numbers  based  upon  the  diameter  between  60,000  and  500,000.  Their 
expression  rewritten  in  terms  of  the  same  notation  as  the  above 
equations  is  given  by 
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Their  use  of  the  jet  to  freestream  mass  flux  ratio  as  opposed  to  the 
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dynamic  pressure  (or  equivalently,  the  momentum  flux)  ratio  provides 
better  agreement  than  many  other  formulae.  Margason  (Ref  9) 
experimentally  determined  another  relationship  describing  the  trajectory 
of  the  circular  air  jet  in  terms  of  the  jet  injection  angle  and 
effective  velocity  ratio.  His  comparisons  with  both  experimental  data 
and  empirical  trajectory  equations  of  other  investigators  provide 
invaluable  insight  into  the  jet  plume  deflection  phenomenon.  Written  in 
terms  of  the  same  variables  as  the  preceding  equations,  Margason's 
relationship  is  expressed  by, 
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where  Vg  is  the  effective  velocity  ratio. 
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is  the  jet  injection  angle. 


Not  only  were  his  comparisons  of  existing  formulae  to  experiment  of 
prime  importance,  but  also  his  determination  through  experiment  that  the 
trajectory  equations  could  be  applied  to  jets  injected  either  normally 
upwards  or  downwards  (for  the  range  of  Vg  considered)  increased  the 
utility  of  the  existing  empirical  data  base.  Nowhere  else  are  data 
ecorded  for  downwardly  injected  jets;  additionally,  his  research  is 
unique  in  that  data  were  presented  for  injection  angles  of  30  through 
180  degrees  (opposing  flow)  and  for  numerous  values  of  effective 
velocity  ratio. 
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Although  the  preceding  trajectory  equations  were  initially 
developed  for  round  jets  injected  into  crossflows,  the  same  expressions 
are  valid  as  first  approximations  for  orifice  shapes  other  than  circular 
if  the  diameter  of  the  jet,  d,  is  replaced  by  the  hydraulic  diameter 
which  is  defined  as  follows: 


cross  sectional  area 


wetted  perimeter 


(10) 


Ivanov,  as  referenced  in  the  text  by  Abromovich,  confirmed  this  fact  for 
rectangular  jets  with  an  aspect  ratio  (length/width)  of  5.0  both  in  the 
streamwise  and  blunt  orientations;  however,  no  reference  was  made  to  the 
flow  conditions  upon  which  this  conclusion  was  based.  It  should  be 
noted  that  a  similar  substitution  occurs  in  convective  heat  transfer 
problems  involving  tubes  of  non-circular  cross  section  (Ref  10) . 

The  effects  of  orifice  shape  on  jet  flow  coefficients  for  normally 
discharging  jets  into  a  cross  flowing  air  stream  were  investigated  by 
Callaghan  and  Bowden  (Ref  11).  This  work  was  expanded  by  Rugger i, 
Callaghan,  and  Bowden  (Ref  12)  to  determine  the  penetration  of  air  jets 
directed  perpendicularly  to  the  free  stream  air  from  non-circular  jets 
and  to  correlate  the  results  with  those  for  circular  jets  through  the 
use  of  the  jet  flow  coefficient.  Ruggeri  (Ref  13)  later  furthered  the 
efforts  of  his  co-workers  by  similarly  correlating  temperature  profiles 
of  heated  air  jets  downstream  of  the  injection  point. 

As  discussed  to  this  point,  each  of  the  investigative  efforts  has 
provided  the  means  to  describe  the  trajectory  of  the  jet-crossflow 
interaction,  but,  none  have  provided  any  information  regarding  the 
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observed  gas  dynamic  phenomena  to  include  the  creation  and  transport  of 
the  counterrotating  vortex  pair  and  the  entrainment  of  both  crossflow 
mass  and  momentum.  Jordinson  (Ref  14)  was  among  the  first  researchers 
to  provide  detailed  data  about  the  structure  of  both  the  flowfield  and 
the  jet  plume.  He  provided  total  pressure  coefficient  contours  and  jet 
centerline  trajectories  for  jet  to  crossflow  velocity  ratios  of  4.  6, 
and  8.  Additionally,  he  presented  a  flowfield  vector  plot  for  3  jet- 
crossflow  velocity  ratio  of  6.  In  his  data,  the  characteristic  kidney 
shape  cross  sections  in  planes  normal  to  the  trajectory  were  clearly 
observed.  The  results  also  show  the  entrainment  of  low  velocity  air 
from  the  vicinity  of  the  injection  plane  into  the  jet.  Gelb  and  Martin 
(Ref  15)  recorded  pressure  distributions  and  conducted  a  detailed 
photgraphic  study  of  high  speed  jets  (125  to  1000  feet  per  second)  being 
injected  into  low  speed  (15  to  60  feet  per  second)  crossflows. 

Keffer  and  Baines  (Ref  16)  formulated  an  entrainment  coefficient 
which  was  proportional  to  the  difference  between  the  jet  and  crossflow 
velocities  and  which  considered  only  the  effect  of  the  entrained 
momentum  upon  the  trajectory  of  an  initially  round  turbulent  non-buoyant 
jet.  Slawson  and  Csanady  (Ref  17)  observed  power  station  chimney  stacks 
and  the  effluent  discharged  from  them.  From  careful  observations  and 
known  values  of  exhaust  temperature,  velocity,  and  mass  flow  in  addition 
to  known  ambient  parameters,  they  formulated  a  modified  turbulent 
entrainment  hypothesis  which  enabled  them  to  predict  the  mean  path  of 
the  turbulent,  buoyant  exhaust  plumes.  Pratte  and  Baines  (Ref  18)  used 
a  methodology  similar  to  Kargason  including  the  use  of  photographs  of  a 
round  turbulent  jet  injected  with  an  aerosol  (oil  as  opposed  to  water) 
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to  obtain  not  only  trajectory  data,  but  also  data  relating  to  the  jet 


spread  and  eddy  structure  of  the  plume.  Chan  and  Kennedy  compared 
numerous  empirical  expressions  (Ref  2,  Table  2-1)  to  compute  the  mean 
trajectory  of  the  jet-crossflow  interaction  for  both  entrainment  and 
drag  formulations,  some  of  which  have  been  discussed  within  this 
section.  They  observed,  as  did  Margason,  that  nearly  all  the  empirical 
equations  defining  the  jet  trajectory  can  be  written  in  the  following 
general  form. 


x 

d 


(11) 


where  K  is  some  prescribed  constant,  and, 

3  and  1  <  n  <  1.5  • 

For  completeness,  Lee  (Ref  19)  and  Garner  (Ref  20)  authored  survey 
reports  detailing  the.  state-of-the-art  of  jet-crossflow  interaction 
studies. 

Experimentalists  of  the  next  decade  carried  on  from  where  those  of 
the  1960's  stopped.  Kamotani  and  Greber  (Ref  21)  not  only  extended  the 
pioneering  work  of  Keffer  and  Baines  for  unheated  circular  air  jets,  but 
they  also  extended  the  available  data  base  for  circular  heated  air  jets. 
The  dependency  of  the  trajectory  on  the  jet  to  crossflow  momentum  flux 
(or,  dynamic  pressure)  ratio,  Q,  was  confirmed.  Additionally,  they  also 
confirmed  that  entrainment  of  the  crossflow  was  independently  controlled 
by  the  normal  and  parallel  components  of  rhe  velocity  impacting  the 
plume,  an  assumption  used  by  Chan  and  Kennedy  as  well  as  other 
investigators-. 
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Experimental  and  analytical  "far  field"  studies  of  viscous  three- 

dimensional  turbulent  air  jets  exhausting  from  rectangular  and  other 

shaped  orifices  into  quiescent  environments  were  conducted  by  Sforza, 

Steiger,  and  Trentacoste  (Ref  22),  Trentacoste  and  Sforza  (Ref  23),  and 

Sforza  Ref  (24).  Detailed  "near  field"  velocity  and  shear  stress 

measurements  were  made  for  turbulent,  incompressible  air  jets  ejected 

from  square  orifices  into  quiescent  environments  by  duPlessis,  Wang,  and 

Kahawita  (Ref  25).  Their  results  showed  square  jet  characteristics 

which  differed  only  slightly  from  those  of  corresponding  circular  jets. 

Only  recently,  have  the  properties  of  non-circular  jets  being  ejected 
» 

normally  into  subsonic  crossflows  been  of  interest.  Weston  and  Thames 
(Ref  26)  conducted  experimental  investigations  of  rectangular, 
isothermal  air  jets  normally  injected  into  subsonic  cross  flowing  air 
from  an  orifice  with  an  aspect  ratio  of  4.  Both  the  streamwise  and 
blunt  orientations  of  the  jet  were  examined  for  both  planar  and 
streamlined  body  (nacelle)  injection  surfaces.  Empirical  trajectory 
power  law  expressions  somewhat  similar  to  those  obtained  for  round  jets 
were  obtained  after  extensive  numerical  analyses  to  ensure  the  best  fit 
for  the  range  of  jet  to  crossflow  velocity  ratios  tested  (4,  8,  and  10). 
Additionally,  empirical  power  law  expressions  similar  to  those  for  the 
trajectory  were  presented  to  define  the  vortex  locus  curve,  i.e.,  the 
trajectory  the  induced  vortices  follow. 

Potential  Flow  Methods.  Many  potential  flow  methods  have  been 
devised  to  ascertain  the  aerodynamic  effects  of  normal  circular  jet 
injection  upon  V/STOL  aircraft.  The  earliest  potential  method  was 
developed  by  Chang  Hsin-Chen  (Ref  27)  to  mathematically  describe  the 
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vortex  "roll-up"  effect  observed  when  a  cylinder  is  placed  in  a  uniform 
crossflow.  Since  jet  injection  from  V/STOL  aircraft  occurs  at 
relatively  high  ratios  of  jet  to  crossflow  velocities,  the  jet  appears 
as  a  quasi-solid  body  and  consequently,  the  methodology  which  Chang 
derived  has  been  used  as  the  basis  for  more  complex  potential  flow 
algorithms.  Bradbury  and  Wood  (Ref  28)  not  only  experimentally  measured 
the  pressure  field  on  the  injection  surface  from  which  a  circular  jet 
issued  for  jet  to  crossflow  velocity  ratios  of  2,  4,  8,  and  11.3,  but, 
they  also  noted  that  the  jet  could  be  qualitatively  simulated  by  a 
blockage  and  distributed  line  sink  representation.  Gelb  and  Martin  (Ref 
15),  attempted  to  describe  the  jet  flow  interaction  for  both  low  and 
high  jet  to  crossflow  velocity  ratios  by  the  superposition  of  a  uniform 
stream  and  a  line  sink  distribution  on  the  jet  axis  with  only  limited 
success.  In  the  same  vein,  albeit  a  more  successful  one,  Wooler  (Ref 
29)  developed  his  own  expression  for  the  jet  plume  centerline  based  upon 
the  experimental  results  of  Jordinson  from  which  the  local  radius  of 
curvature  along  the  jet  axis  could  be  determined.  After  balancing  the 
centrifugal  force  exerted  on  the  jet  by  the  crossflow,  the  jet  plume  was 
subdivided  into  a  large  but  finite  number  of  elements.  A  distribution 
of  trailing  vortices  was  then  developed  from  which  the  induced  flow 
field  and  pressure  distribution  on  the  injection  surface  could  be 
computed.  Williams  and  Wood  (Ref  30)  developed  a  more  refined  "vortex 
sheet"  method  based  upon  the  earlier  work  of  Wooler.  It  also  required 
an  empirical  expression  for  the  jet  trajectory,  but,  it  provided  better 
agreement  with  experimental  data.  Wooler,  Burghart,  and  Gallagher  (Ref 
31)  developed  a  theoretical  model  which  did  not  require  the  use  of  an 


13 


empirical  expression  for  the  jet  center  line,  but  which  required  the  jet 
to  be  represented  by  a  distribution  of  both  sinks  and  doublets  to 
account  for  the  entrainment  and  blockage  of  the  jet.  Again,  the 
objective  of  the  investigators  was  to  determine  the  induced  aerodynamic 
effects  upon  the  surface  from  which  the  jet  was  issuing  and  not  to 
determine  quantitative  data  concerning  the  physics  of  the  flowfield 
interaction.  Wu,  Mosher,  and  Wright  (Ref  32)  conducted  both  analytical 
and  experimental  investigations  for  circular  and  elliptical  normally 

injected  jets  in  streamwise  and  blunt  configurations.  Their  potential 
method  utilized  doublets,  vortices,  and  sinks  to  model  the  effects  of 
entrainment  and  blockage;  additionally,  the  wake  region  of  the  jet  was 
excluded.  However,  this  model,  as  well  as  other  potential  flow  models, 
required  the  use  of  an  empirical  formula  for  the  jet  trajectory  to 
describe  the  three-dimensionality  of  the  flow  interaction  which  simply 
is  not  available  in  potential  representations.  The  empirical  model 
utilized  in  their  method  was  that  of  Wooler  which  is  presented  here  for 
completeness, 

x  /  y  v 

-  -  B  ( Cosh  -  -  1  )  (12) 

d  V  Bd  ' 

where  B  ■  0.19 

Additional  investigations  have  followed  similar  approaches.  Among  them 
are  the  works  of  Braun  and  McAllister  (Ref  33)  and  Rubbert  (Ref  34). 
Endo  and  Nakamura  (Ref  35)  developed  a  mathematical  model  of  the  round 
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jet-crossflow  interaction  by  consolidating  the  efforts  of  previous 
investigators,  including  the  vortex  model  of  Wooler,  and  then  correlated 
predictions  from  it  with  experimental  data. 

Integral  Methods.  The  work  of  Chan  and  Kennedy  included 
approximate  integral  analyses  of  the  turbulent  round  jet  for  all  three 
zones  of  the  jet-crossflow  interaction  as  well  as  experimentally  derived 
data.  Comparisons  of  their  analytical  predictions  with  experimental 
data,  both  their  own  as  well  as  others,  consistently  show  reasonable 
agreement.  Their  bibliography  of  jet-crossflow  research  efforts  is  by 
far  one  of  the  most  complete  and  current.  Schatzmann  (Ref  36)  extended 
the  efforts  of  Chan  and  Kennedy  for  turbulent  round  jets  in  both  the 
"curvilinear"  and  the  "far  field"  regions  and  corrected  the  mathematical 
errors  entailed  in  the  analyses  from  which  Chan  and  Kennedy  derived 
their  methodology.  Additionally,  Schatzmann  recast  the  governing 
equations  in  terms  of  the  natural  coordinates  of  the  jet  in  vector  form. 

Numerical  Methods.  Numerical  computation  of  the  three-dimensional 
behavior  of  jets  has  become  possible  with  the  current  generation  of  high 
speed  computers.  Computations  of  co-flowing  turbulent  supersonic 
circular,  elliptic,  and  interacting  square  jets  were  accomplished  by  Oh 
and  Harris  (Ref  37)  utilizing  an  ADI  (Alternating  direction  implicit) 
scheme  to  solve  the  steady  flow  parabolic-elliptic  Navier-Stokes 
equations  in  three  dimensions.  The  parabolized  equations  were  formed 
from  the  Navier-Stokes  equations  utilizing  the  assumption  of  a  preferred 
main  flow  direction  which  dominates  the  convective  processes. 
Consequently,  derivatives  with  respect  to  this  direction  are  assumed 
negligible  in  the  shear  stress  terms.  However,  as  noted  by  Oh  and 
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Harris,  "...the  system  is  formally  correct  for  supersonic  flows;  for 
subsonic  flows,  the  static  pressure  distribution  or  its  gradient  in  the 
marching  direction  must  be  specified  if  the  problem  is  to  be  well  posed 
mathematically".  Little  additional  complexity  was  added  to  their 
methodology  by  the  incorporation  of  a  simple  algebraically  prescribed 
mixing  length  eddy  viscosity  model  through  which  turbulence  closure  was 
achieved.  McGuirk  and  Rodi  (Ref  38)  computed  the  steady,  three- 
dimensional  behavior  of  turbulent,  rectangular,  incompressible,  free 
jets  with  aspect  ratios  of  1,  5,  10,  and  20  issuing  into  a  quiescent 
atmosphere.  They  utilized  the  finite  difference  scheme  of  Pantankar  and 
Spalding  for  three-dimensional  parabolic  flows  and  achieved  closure 
through  use  of  a  two-equation  turbulent  kinetic  energy-dissipation  (k-€) 
model.  Chien  and  Schetz  (Ref  39)  computed  the  normally  injected 
turbulent  round  jet-crossflow  interaction  for  an  incompressible  fluid 
using  the  vorticity-velocity  formulation  of  the  steady  Navier-Stokes 
equations.  Closure  was  obtained  through  use  of  a  Prandtl  constant  eddy 
viscosity  model.  Pantankar,  Basu,  and  Alpay  (Ref  40)  applied  a  finite 
difference  code  to  solve  the  turbulent  round  jet  injection-crossflow 
problem  for  an  incompressible  fluid.  Closure  was  achieved  through 
essentially  the  same  two-equation  turbulent  kinetic  energy-dissipation 
model  McGuirk  and  Rodi  used.  The  steady,  incompressible  formulations  of 
the  last  two  groups  of  investigators  numerically  predicted  the  flowfield 
caused  by  the  interaction  of  the  jet  and  crossflow  for  relatively  slow 
crossflows  and  high  jet  to  crossflow  velocity  ratios  with  reasonable 
accuracy. 
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1.2  Summary 

Numerous  empirical  formulae  have  been  developed  to  describe  the 
trajectory  of  a  jet  injected  at  some  arbitrary  angle  into  a  crossflow 
comprised  of  the  same  or  a  different  fluid.  Most  of  these  expressions 
are  generalizable  to  a  power  law  form  as  expressed  by  Eq  (11);  however, 
these  expressions  were  developed  for  flows  characterized  by  neither  low 
ratios  of  jet  to  freestream  momentum  flux,  nor,  for  low  ratios  of  jet  to 
freestream  mean  molecular  weights.  Potential  flow  approximations  to  the 
jet-crossflow  interaction  problem  have  marginal  utility  for  flowfield 
analysis  although  they  may  be  used  to  adequately  describe  the 
aerodynamic  effects  of  the  jet  upon  V/STOL  aircraft  performance, 
stability,  and  control.  Numerical  methods  currently  being  used  to  solve 
the  three-dimensional  jet  injection  problem  are  restricted  to  jets 
injected  into  relatively  slow  crossflows  or  quiescent  atmospheres. 
Furthermore,  the  Navier-Stokes  equations  were  simplified  by  assuming  no 
time  dependence  and  constant  fluid  density.  In  some  cases,  the  set  of 
equations  was  parabolized  to  enable  existing  solution  algorithms  to  be 
used.  To  date,  no  other  computational  methodology  has  been  documented 
which  will  permit  computation  of  the  the  "omplex  flow  structure  caused 
by  the  interaction  of  real,  viscous  fluids  being  injected  into 
crossflows  characterized  by  either  low  or  high  jet  to  crossflow  velocity 
ratios.  Furthermore,  no  known  method  developed  to  date  has  been  used  to 
even  attempt  the  same  feat  with  real,  low  molecular  weight  gas  mixtures 
at  elevated  temperatures  injected  into  a  high-subsonic  ambient  crossflow 
with  atmospheric  properties  corresponding  to  the  altitude  at  which  the 
injection  occurs.  The  balance  of  this  document  details  the  development 
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and  validation  of  the  methodology  as  a  useful  predictive  engineering 
tool  to  advance  both  the  state-of-the-art  and  the  general  knowledge  of 
jet-crossflow  interactions.  Section  II  describes  the  governing 
equations,  boundary  conditions,  and  ancillary  relations  required  to 
mathematically  model  the  physical  problem.  Section  III  details  both  the 
finite  difference  formulations  of  the  model  and  the  methodology  to 
implement  them  as  well  as  the  factors  inherent  in  both  the  model  and 
physical  problem  which  directly  affect  the  stability  and  convergence  of 
the  numerical  solutions.  Finally,  the  results  of  the  research  and  the 
associated  conclusions  drawn  from  them  are  presented  in  Sections  IV  and 
V  respectively. 


II  *  Mathematical  Model  of  the  Physical  Problem 

The  complex  flow  interaction  problem  created  by  the  exhaust  of  an 
airborne  chemical  laser  aboard  a  wide-bodied  aircraft  which  is 
discharging  into  a  high-subsonic  ambient  crossflow  presents  formidable 
difficulties  in  the  modeling  of  the  mixing  processes.  The 
simplifications  made  by  previous  investigators  to  permit  integral 
analyses  to  be  conducted  are  not  valid  for  the  existing  problem  and, 
consequently,  the  complete  set  of  non-lLnear,  coupled  equations  (Navier- 
Stokes)  must  be  utilized  to  describe  the  behavior  of  the  visccus- 
inviscid  fluid  interaction;  the  equations  are  further  encumbered  by  the 
requirement  to  compute  the  fluid  properties  of  the  mixture  which  are 
constantly  changing  with  the  temperature  and  species  concentration. 
These  equations  are  formulated  in  Cartesian  tensor  notation  for 
compactness  and  clarity  within  this  section. 

Additionally,  the  wide  variation  in  molecular  weights  further 
exacerbates  the  task  of  approximating  properties  of  the  fluid  mixture. 
Using  the  assumption  that  the  exhausting  multi-component  gas  is 
homogeneous  in  composition  while  neglecting  the  possibility  of  chemical 
reactions  with  ambient  air  occurring  within  the  fioafield  permits  not 
only  the  binary  mixture  diffusion  coefficient  (Appendix  C)  to  be 
computed,  but,  also  introduces  many  simplifications  in  both  the  species 
conservation  and  global  energy  conservation  equations. 

Because  the  matrix  of  problem  dependent  parameters  is  overwhelming 
in  nature,  representative  values  for  the  aircraft  Mach  number  (0.7)  and 
pressure  altitude  (36089.239  feet,  or,  11  kilometers)  were  chosen. 
Additionally,  only  one  representative  jet  mixture  (Table  I)  was 
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utilized.  For  all  real  gas  test  cases,  the  jet  temperature  was  defined 
to  be  1500  °K  and  the  jet  Mach  number  was  defined  to  be  0.24;  the  ratio 
of  jet  to  freestream  dynamic  pressures  was  constrained  to  the  value 
0.15.  The  orifice  aspect  ratio  of  1-.75  used  for  the  injection  studies 
was  based  on  preliminary  design  data  of  the  laser  exhaust  diffuser  as 
implemented  on  the  bottom  of  the  fuselage  of  a  wide-bodied  air  vehicle. 


II. 1  Model  Governing  Equations 


The  full  set  of  time  dependent  Navier-Stokes  equations  describing 
the  mean  flow  quantities  in  three  dimensions  in  addition  to  a  species 
conservation  equation  were  required  to  solve  the  jet  injection  and 
associated  flowfield  interaction  problem.  These  equations  for  the  gas 
mixture,  in  conservation  form,  are  expressed  in  Cartesian  tensor 
notation  (utilizing  the  Einstein  summation  convention)  as  follows: 


Conservation  of  Mass 


P»j-  +  (P  Uj  ) » j 


Conservation  of  Momentum 


(P  «i  ),t  +  (  PU^  +  P,±  -  -  pfi 


»h«.  Ty  -X  «  «  +  Vef£«  °1>3  +  U3,l  > 


and  where  ^  «  -  -  U 


-  (  V  +wt  ) 
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Conservation  of  Energy 


(  PE  ),t  +  (  P EUj  +  PUj  -  Tkjuk  +  qj  ),j  -  Pf^  !38) 

where  E  -  I  +  {  U±  V±  )/2  (19) 

N 

and  where  qj  «  -  keffTtj  +^Gdiff  k  hk  (20) 

k 

where  6dlff  k  -  -  P&eff  (  Pk/  P  ),j  (2.'-a) 

keff  -  k  +  »Jt  Cp/Prt  ^21“b> 

Peff  “  V  +  (  *VP)/Sct  l21-c) 


Conservation  of  Species 

<Pk  ),t  +  (PfcOj  +  Gdlff  k  ).J  -  0  (22) 

Despite  the  favorable  properties  and  utility  of  the  conservation  form  of 
the  governing  equations,  sometimes  if.  is  even  more  desirable,  depending 
upon  the  solution  algorithm,  to  reformulate  them  in  non-cons ervati'”  , 
primitive  variable  form.  This  was  accomplished  on  the  momentum 
equations  by  expanding  Eq  (14)  and  then  removing  the  continuity  equation 
(  Eq  (13)  ).  The  resultant  expression  for  the  momentum  equations  is 
given  by. 


p(Vt  +  Vifj  >  +  p’i  -  lij.j  -  pfi 


(23) 
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Similarly,  the  continuity  equation  was  removed  from  the  energy  equation. 
Specifically,  when  Eq  (23)  is  dotted  with  the  velocity,  U^,  the 
following  expression  results: 


Pi  (  UjUi/2  ),t  +  11^  U^u/2  }ti  ]  + 


UiP,i  “  UiTi.1-j  “  pQifi  *  0 


This  expression  was  then  subtracted  out  of  the  energy  equation  to 
provide  the  following  ancillary  version  of  the  energy  equation  which  was 
acrually  solved. 


P(  I.t  t  OjX.j  )  +  PU^J  -  Tkj0k(.  *  htl  -  0  . 


Thu:.,  the  Navier-Scohes  equations,  viz.  das  (13),  (23)  and  (25),  and  the 
species  conservation  equation,  Eq  (22),  comprise  the  set  of  governing 
equations  of  the  model. 


L'L . 2  Therraodynf.tuic  Properties 

The  thermodynamic  properties  of  the  gas  mixture  are  calculated  as 
functions  of  temperature  frov-  tabular  data  of  each  constituent  gas 
(Appendix  2)»  Both  the  specific  heat  at  constant  p*srersure,  Cp,  and 
the  enthalpy  of  the  gas  mixture,  h,  are  calculated  from  the  following 
equations : 


CP  Zj  <pj;p  > c?,  ■  Ebj  CI 

j  ■  j 

*■  ■  E  »j » j 
j 
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where  the  expression  (  P^/P  )  is  the  mass  fraction 


or  tne 


j  -species,  m. ,  and  where  and  h  r.re  tha 

j  Pj  J 

jtu-species  specific  heat  at  constant  pressure  end 
enthalpy,  respectively. 


.th 


Tile  average  molecular  weight  of  the  gas  mixture,  M,  i3  defined 
(Ref  10)  by  the  following  expression, 


(28) 


where  is  the  j  ^-species  molecular  weight.  The  ratio  of  specific 
heats  at  constant  pressure  to  that  at  constant  volume,  Y  ,  is  then 
expressible  in  terms  of  the  gar  mixture  variables  which  were  presented 
above,  that  is. 


Y  =  cp  /  cv  -  jl  -  Rgas  /  H  Cp  |  (29) 

The  equation  of  state  for  an  ideal  gas  was  written  in  terms  of  the 
specific  internal  energy,  I,  and  the  ratio  of  specific  heats  (Eq  (29)); 
it  is  given  by  the  following  expression! 

P  -  (  Y  -  1  )  PI  (30) 

The  temperature  of  the  gas  mixture  was  then  calculated  locally  by  the 
following  expression. 


1  1  Sr 


iy  / 


(31) 
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Finally,  the  local  speed  of  sound  of  the  gas  mixture  was  calculated  from 
the  following  equation: 


Y(  Y  -  1  )  I 


II. 3  Turbulence  Model 

To  incorporate  the  effects  of  turbulence,  a  simple  algebraically 
prescribed  eddy  viscosity  model  was  used  instead  of  any  higher  order 
model  which  did  not  necessarily  improve  the  accuracy  of  the  solution 
while  it  most  assuredly  would  increase  the  computational  complexity. 
Harsha  (Ref  41)  reviewed  numerous  locally  and  history  dependent  jet 
turbulence  models  for  widely  varying  gas  combinations  in  co-flowing 
jets;  however,  none  were  evaluated  explicity  for  jets  discharging 
normally  into  crossflows.  Chien  and  Schetz  (Ref  39)  used  a  simple 
constant  Prandtl  velocity  defect  eddy  viscosity  formulation,  the  value 
of  which  was  predetermined  by  the  initial  jet  and  crossflow  velocities. 
They  further  determined  that  for  an  incompressible  fluid,  variations  in 
the  proportionality  constant  in  this  model  affected  the  jet  trajectories 
minimally  and  therefore  retained  the  constant  associated  with  coflowing 
round  jets  such  that  the  eddy  viscosity  was  expressed  by  (Ref  39,  Eq  16) 


0.0256 


(  Ujet  -  Va  )  (  D/2  ) 


where  D  is  the  injection  diameter  of  the  jet.  Morkovin  (Ref  42)  has 
shown  that  the  turbulence  structure  of  gasdynamic  flows  is  essentially 
unaffected  by  compressibility  for  Mach  numbers  less  than  5.  Because  of 


the  complexity  of  the  flowfield  and  the  inherent  difficulty  in  modeling 
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it,  the  simplest  model  to  reasonably  approximate  the  turbulent  transport 
phenomena  should  be  used.  Consequently,  Eq  (33)  was  modified  to  account 
for  variations  of  jet  and  freestream  velocities  as  well  as  mixture 
density.  The  resultant  expression  for  the  turbulence  nodal  used  in  the 
flowfield  simulation  is  given  by, 

y_  -  0.0256  P±  |U  -  UJ  <  D/2  )  (a4) 

where  U  «  y/  Ui  Ui  ,  i  *  i,  2,  3 

and  is  the  local  mixture  density. 

Additionally,  Eq  (34)  has  the  required  property  of  applying  the 
turbulence  corrections  where  the  jet  i3  actually  mixing  with  the 
crossflow.  As  the  flow  becomes  more  uniform,  the  augmented  mixing  due 
to  turbulence  becomes  negligible  and  the  mean  flow  (laminar)  equations 
correctly  define  the  fluid  dynamics. 

II. 4  Boundary  Conditions 

The  approximation  of  the  physical  problem  was  accomplished  by 
creating  a  rectangular  computational  domain  (Figure  3).  The  faces  of 
this  domain  are  defined  as  follows: 

(1) .  Upper  face 

a.  Injection  surface,  "no-slip"  boundary 

b.  Jet  orifice,  prescribed  inflow  boundary 

(2) .  Windward  face,  prescribed  inflow  boundary 

(3) .  Leeward  face,  continuative  outflow  boundary 

(4) .  Front  lateral  face,  continuative  outflow  boundary 

(5) .  Back  lateral  face,  continuative  outflow  boundary 
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(6).  Lower  face,  confirmative  outflow  boundary 

The  rectangular  jet  was  forced  to  discharge  into  a  uniform  crossflow 
with  ambient  properties  corresponding  to  the  aircraft  pressure  altitude; 
the  local  boundary  layer  build-up  corresponded  to  that  of  a  wide  bodied 
aircraft  at  the  proposed  location  for  the  chemical  laser  diffuser  exit* 
The  capability  to  prescribe  either  ’'top-hat"  or  Gaussian  (Appendix  B) 
profiles  of  the  jet  injection  velocity  was  accommodated  in  the  model. 
The  boundary  conditions  employed  in  the  mathematical  model  are  presented 
as  follows  for  each  face  of  the  computational  domain. 

II. 4.1  Upper  Face  (6) 

The  injection  surface  of  the  upper  face  (simulating  the  lower 
surface  of  the  aircraft  fuselage)  was  defined  to  be  "no-slip"  such  that 
both  the  normal  and  tangential  components  of  the  velocity  were  zero, 
that  is. 


Ui  -  0  ,  i  =  1,  2,  3  .  (35) 

It  was  further  assumed  that  the  injection  surface  behaved  as  an 
adiabatic  wall.  This  condition  results  in  derivatives  of  the  fluid 
specific  internal  energy,  I,  at  the  surface  in  the  normal  direction 
exhibiting  no  variation  assuming  that  cy  is  locally  constant  for  every 
affected  cell.  When  the  boundary  layer  approximation  of  the  pressure 
variation  in  the  normal  direction  is  also  incorporated  within  the  model, 
the  following  general  fluid  property  boundary  condition  results  at  the 
injection  surface,  namely. 
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0 


(36) 


where  P  is  any  thermodynamic  variable  of  the  fluid, 
such  as  pressure,  temperature,  internal  energy,  etc. 


The  prescribed  inflow  conditions  of  the  jet  orifice  were  defined  by  the 
properties  of  the  exhaust  gas  mixture  (Table  I)  at  an  arbitrarily 
prescribed  temperature  (1500  °K)  and  Mach  number  (0.24).  The  jet  to 
crossflow  dynamic  pressure  ratio  was  constrained  to  be  0.15  which  in 
turn  specified  the  exhaust  jet  pressure  and  density  since  both  the 
freestream  and  jet  Mach  numbers  were  held  fixed. 


II. 4. 2  Windward  Face  (1) 

This  face  of  the  computational  domain  was  a  completely  prescribed 
inflow  boundary  whose  fluid  properties  and  velocities  were  defined  by 
the  aircraft  flight  parameters  of  Mach  number  and  altitude.  The  normal 
(streamwise)  component  of  velocity  was  prescribed  by. 


Uoo(y/^)  »  within  the  turbulent  boundary 

layer,  where  the  boundary  layer 
thickness  is  given  by. 


0.376  x  Re. 


,  outside  the  turbulent  boundary  layer. 


whereas  the  tangential  velocity  components  were  both  assumed  to  be 
negligible,  that  is. 


v  *  w  *  0  . 
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Again,  the  normal  derivative  of  the  fluid  thermodynamic  variables  was 
assumed  to  be  zero,  and  specified  mathematically  by  Eq  (36). 

11.4*3  Remaining  Faces  (2 ,3, 4, 5) 

The  remaining  faces  of  the  computational  domain  were  all  defined  as 
continuative  outflow  boundaries.  Initially,  the  normal  derivatives  of 
both  the  velocity  and  the  gas  mixture  properties  were  enforced  to  be 
identically  zero,  i.e., 

(9U.  )  -  0  (39) 

9n 

and 

(  9 P  )  -  0  ,  (36) 

9n 

where  P  again  is  any  fluid  thermodynamic  variable.  However,  as  the 
solution  of  the  equations  progressed,  the  velocities  and  properties  were 
allowed  to  adjust  to  ensure  compliance  with  both  the  boundary  conditions 
and  the  governing  equations. 

II. 5  Summary 

The  physical  problem  of  an  airborne  chemical  laser  venting  waste 
exhaust  gases  overboard  into  a  high-subsonic  (0.7  Mach)  ambient 
crossflow  at  some  prescribed  aircraft  pressure  altitude  (11  km)  was 
modeled  by  creating  a  rectangular  computational  domain  over  which 
certain  necessary  boundary  conditions  were  enforced.  The  upper  surface 
of  this  domain  was  considered  to  represent  the  "no-slip"  aircraft 
fuaelage  bottom  through  which  the  rectangular  exhaust  jet  was  forcibly 
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ejected  downwards  at  a  given  temperature  (1500  °K),  Mach  number  (0.24), 
and  jet  to  crossflow  dynamic  pressure  ratio  (0.15).  The  remaining  faces 
were  established  as  continuative  outflow  boundaries.  Only  the  leeward 
face,  normal  to  the  incoming  crossflow,  was  prescribed  by  the  aircraft 
Mach  number  and  pressure  altitude;  the  others  were  permitted  to  relax 
from  the  initial  prescription  of  the  "no-change",  or  "first  derivative", 
boundary  condition  and  to  acquire  values  which  ensured  that  the 
governing  Navier-Stokes  equations  Here  satisfied.  The  next  section 
details  the  finite  difference  algorithm  and  the  requisite  subdivision  of 
the  computational  domain  employed  to  solve  the  mathematical  description 
of  the  model  presented  in  this  section. 


Ill*  Numerical  Solution  Algorithm 


m. 

G$r 


i 


m. 


The  computational  domain  defined  in  Section  II  was  divided  into  a 
large,  but  finite,  number  of  three-dimensional  rectangular  cells  as 
shown  in  Figure  4.  The  governing  non-linear  differential  equations  were 
approximated  by  finite  difference  algebraic  equations  which  were  then 
solved  for  every  cell  within  the  computational  domain.  A  simple 
algorithm  was  developed  by  Hirt,  Nichols,  and  Romero  (Ref  43)  for  time- 
dependent,  two-dimensional,  incompressible  flows  and  was  later  modified 
by  Cloutman,  Nichols,  and  Romero  (Ref  44)  to  incorporate  the  effects  of 
compressibility  with  the  assumption  of  constant  fluid  transport 
properties.  The  latter  version  also  incorporated  the  implicit 
Continuous  fluid  Eulerian  (ICE)  method  as  an  option  (Ref  45)  by  which 
the  best  estimate  of  the  time  advanced  pressure  was  utilized  in  the 
momentum  equations.  The  algorithm  used  a  primitive  variable,  donor-cell 
formulation  which,  in  fact,  formed  the  basis  of  the  methodology  used  for 
the  numerical  simulation  of  the  jet  injection  problem.  However,  the 
algorithm  was  further  modified  in  this  study  to  solve  fully  three- 
dimensional,  mixed  species,  real  gas  flows  with  both  variable  transport 
and  thermodynamic  properties. 

III.l  Computational  Grid 

The  spatial  increments  in  the  x  and  y  directions  respectively,  Ax 
and  '  y,  were  uniform,  that  is, 

I  A x  -  constant, 

J  A-  -  constant,  (40) 

v 
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but,  they  were  not  necessarily  equal.  Because  of  the  violent  nature  of 
the  jet  injection  problem  and  the  associated  near  field  mixing,  the 
corresponding  spatial  increments  in  the  z-direction  were  varied 
geometrically  from  the  injection  surface  in  an  effort  to  capture  the 
physics  of  the  problem  by  locally  refining  the  mesh  in  the  region  of 
maximum  fluid  interaction.  The  mathematical  description  of  the 
geometrically  varying  cell  lengths  in  the  z-direction  is  given  by, 

Azk  -  Azq  Kk_1  (41) 


where  the  AzQ  and  A  are  the  initially  prescribed  and  subsequent  k-th 
values  respectively  of  the  cell  vertical  dimensions  as  measured  from  the 
wall,  and  K  is  the  constant  multiplicative  factor  of  the  geometric 
expansion  between  any  two  consecutive  members. 

Together  with  the  spatial  incrementation  given  in  Eqs  (40)  and 
(41),  the  solution  algorithm  used  a  staggered  mesh  as  indicated  in 
Figure  5.  The  properties  of  the  gas  mixture  for  every  (i, j ,k)-cell,  for 
example,  P(i»j,k),  were  defined  at  the  cell  center;  the  velocities,  on 
the  other  hand,  were  defined  at  the  centers  of  the  (i,j,k)-cell  faces. 
The  velocities  on  the  six  faces  (  u(i,j,k),  u(i-l,j,k),  v(i,j,k),  v(i,j- 
l»k),  w(i,j,k),  w(i,j,k-l)  )  are  indicated  in  the  figure  noting  that  the 
velocity  components  are  situated  in  planes  normal  to  their  respective 
directions . 


III. 2  Spatial  Derivatives  :  Non-uniform  Mesh 


Because  of  the  prescribed  variation  of  A  zk„  the  normal  expressions 


for  the  three-point,  second  order,  spatially  centered  approximations  to 
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the  first  and  second  derivatives  associated  with  a  constant  mesh  could 


not  be  used-  However,  Taylor  series  expansions  were  utilized  to  derive 
the  required  derivative  expressions  for  non-uniform  meshes,  the 
discussion  of  which  follows  directly. 

If  the  set  of  values  of  a  given  parameter  {<{>  ^  >  are  distributed  on 
a  non-uniform  mesh  in  an  arbitrary  spatial  direction  at  the  nodes 
{  ^l»  <**2»  *  *  *  ^n*  *  *  *  \  c^en  the  Taylor  series  expansions  for  both  Vi 
and  may  be  generated  in  terms  of  A£*”  and  AC  ”,  where 


(42) 


4»n-l  -  4>n  -  (AD  +  (AD2 
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-  (AD3  +  ^  (AD4  -••• 

3?  47 


(44) 


When  Eq  (43)  is  multiplied  by  (AD2  and  Eq  (44)  is  multiplied  by 


(A?+)^  ,  and  when  the  resulting  expressions  are  subtracted  from  one 
another,  the  required  expression  for  the  first  derivative  of  <p  is  easily 
obtained  and  it  is  given  by, 

\  **  (AC“)  4>n+1  +  (AC4-  -  AS")  4>n 

to?)  (AC+  +  AO  to?)  (AD 

to?)  ^  +  0(A?+,AD  (45) 

toD  to?  +  AD 

Verification  of  this  spatially  "centered"  expression  for  a  uniform  mesh 
is  easily  obtained  by  requiring  that  A£+  =  AC*  =  A.  The  resulting 
expression, 

VUlti  -  <><a2>  iw 

2A 

is  immediately  discerned  as  the  second  order  accurate,  spatially 
centered  finite  difference  approximation  to  the  first  derivative  of 
for  a  uniformly  varying  computational  mesh. 

Similarly,  the  equivalent  spatially  centered  expression  for  the 
finite  difference  approximation  for  the  second  derivative  of  $  is 
obtained.  When  Eq  (43)  is  multiplied  by  (AC”)  and  Eq  (44)  is  multiplied 
by  (A£+)»  and  the  resulting  expressions  are  added  to  one  another,  the 
required  expression  for  the  second  derivative  of  <J>  is  obtained.  It  is 
given  by. 
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2  *n  + 

(a?TTaO 

Vl  +  OtA^.AO  {47) 

Again,  verification  of  the  expression  for  the  non-uniform  mesh  is 
obtained  by  letting  the  non-uniform  increments  be  equal  to  some 
arbitrary  constant,  that  is,  A£*  «  AC  -A  .  The  resulting  expression, 

C  -  ^n-fl  ~  2  +  +  0(  A2)  <48) 

A2 

is  readily  discerned  to  be  the  spatially  centered,  second  order  accurate 
finite  difference  approximation  for  the  second  derivative  of  <J>  for  a 
constant  mesh  increment.  Both  Eq  (45)  and  Eq  (47)  show  that  typically 
second  order  accurate  finite  difference  approximations  associated  with 
uniform  grids  degrade  drastically  to  first  order  approximations  when  the 
mesh  is  appreciably  non-uniform.  The  larger  the  value  of  the  geometric 
expansion  factor  (K)  in  Eq  (41),  the  more  pronounced  the  numerical  error 
becomes . 

III. 3  Optimum  Cross-Derivative  Approximations 

In  the  momentum  partial  differential  equations,  Eq  (14),  the 
following  velocity  cross  derivative  terms  appear  in  the  expansion  of  the 
stress  tensor: 

x-direction  :  3^  ,  3^ 

3x  3y  3x  3z 


C  - 


n+1 


(h?)  (hC  +  AO 


(AD  (AC*  +  AD 
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y-direction  : 


*mK*pt>*)f\)  m  iffwwwaiwwi 


32u  ,  3^w 


z-direction  : 


(49) 


3y  3x  3y  3z 
32u  9  3^v 


3z  3 


z  ox 


3z  3y 


Each  of  these  partial  derivatives  can  be  approximated  by  suitable  finite 

difference  expressions.  However,  proper  use  of  Taylor  series  expansions 

for  the  velocities  about  the  (i,j,k)-cell  component  velocities  results 

in  approximations  with  negligible  truncation  error. 

2 

As  an  example,  the  approximation  for  9  u 


3z  3x 

in  the  z-direction  momentum  equation  is  fabricated  from  several  Taylor 
series  expansions  about  u(i,j,k).  Since  the  marching  direction  for  the 
expansions  is  identified  by  the  particular  momentum  equation  (z- 
direction) ,  the  u-componenet  in  the  z-direction  is  allowed  to  vary 
between  the  "k"  and  "k+l:i  cells;  because  x  is  not  the  marching 
direction,  the  u-component  is  only  permitted  to  vary  between  the  "i"  and 
"i-1"  cells;  the  y-direction  cell  index,  "j",  remains  fixed. 
Consequently,  the  u-component  velocities  of  interest  are  u(i,j,k), 
u(i,j»k+l),  u(i-l»j ,k+l) ,  and  u(i-l,j,k).  Expansions  of  these 
velocities  about  u(i,j,k)  are  as  follows: 


u(i,j,k+l)  »  u 


+  Az  9u 
i ,  j  ,k  9  z 


+  Az2  92u 

TT  9Z2 


i.j.k 

+  0(Az3) 


(50) 


1»  j  »k 
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,k-fl)  »  u 


+  A?.  3u  j 


•'LX  3u  | 


i,j,k 


3z|;Lj,k  ,  gxli.j.k 


+  ;( 


Ax2  32u 


ox 


-  2  Ax  Az  S2a  +Az232u 


3z  3x 


d  l»  1  \ 

)  z"Ji,i,Kr 


+  0(-Ax3,Az3)  (51) 


u(i-l,.i,k)  «  aj 

-Ax  3u  | 

+  Ax^  3  ^u 

j  >k 

3x  1 

i,j,K  ^x2 

-  0(Ax3) 


(52) 


When  Eqs  (50)  and  (52)  are  added,  and  then  Eq  (51)  is  subtracted  from 
the  resulting  expression,  the  "optimum"  finite  difference  form  of  the 
partial  derivative  is  given  by 


o2u  -  u(i,j,k+l)  -  u(i“l,j,k+l)  -  u(i,  j  ,k)  +u(i-l,j,k.)  (53) 
3z  3x  A  z  A  x 

This  form  of  the  finite  difference  approximation  to  the  example  partial 
derivative  possesses  the  fortuitous  property  of  truncation  error 
cancellation  by  the  components  thereby  increasing  the  accuracy  of  the 
numerical  solution.  The  remaining  finite  difference  approximations  to 
the  partial  derivatives  Indicated  in  Eq  (49)  were  similarly  determined 
and  are  summarized  as  follows: 


x-direction  ; 


32v  -  v(i+l,j,k)  -  v(i+I,j-l,k)  -  v(i> j >k)  +  v(i,j-l,k)  (54) 

9x  9y  A  xA  y 


9*w  «  w(l*kl , j ,k)  -  w(iMtj,k-l)  -w(i,j,k)  +w(i,j,k-l)  (55) 

9*  9s  A  x  A  z 
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y-direction  :  , 


?2u  -  u(L,j+l,k)  -  u(i-l,j+l,k)  -u(i,j,k)  +u(i-l>j,k)  (56) 

8y  8x  Ay  Ax 

a2w  -  w(i,j+l,k)  -  w(i,j+l,k-l)  -wd.j^)  -fw(i,j,k-l)  (57) 

3y  3z  Ay  A3 

z-direction  : 


92v  =  v(i,j,k+l)  -  v(i,j-l,k+l)  -  v(i,j,k)  -  v(i,j-l,k)  (58) 

3z  9y  Az  Ay 

III. 4  Finite  Difference  Equations 

In  accordance  with  the  SOLA-ICE  algorithm  (Ref  44),  the  governing 
partial  differential  equations  were  approximated  by  finite  difference 
equations  using  a  modified  cicnor  cell  formulation.  The  advected  terms 
were  modeled  by  a  sum  of  spatially  centered  and  donor  cell  terms,  while 
the  remaining  terms  were  approximated  by  difference  expressions  and 
techniques  readily  available  in  the  literature.  The  algorithm  utilized 
the  concept  of  a  staggered  mesh  in  which  the  velocities  are  cell  face 
centered  and  the  properties  of  the  fluid  mixture  are  cell  centered. 
This  formulation  permits  the  same  computational  index  to  be  employed  for 
both  the  cell  properties  and  the  cell  component  velocities.  However, 
particular  care  must  be  exercised  during  the  construction  of  the 
difference  equations  because  the  spatial  locations  of  the  computed 
variables  within  each  cell  are  not  coincident  (see  Figure  5)-  When  the 
mesh  incrementation  is  permitted  to  become  infinitesimally  small,  the 
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finite  difference  approximations  to  the  governing  partial  differential 
equations  are  consistent. 

In  FORTRAN  IV,  fractional  and  zero  indices  are  prohibited. 
Consequently,  the  indices  for  the  cell  velocities  were  redefined.  If 
the  (i,j,k)-cell  is  specified,  then  the  properties  of  the  fluid  are 
associated  with  the  spatial  coordinates  of  the  cell  center;  the  x- 
direction  velocity  is  specified  at  the  center  of  the  right  cellular 
face;  the  y-direction  velocity  is  specified  at  the  center  of  the  back 
cellular  face;  the  z-direction  velocity  is  specified  at  the  center  of 
the  upper  cellular  face.  The  velocities  at  the  centers  of  the  remaining 
cellular  faces  are  identified  by  decrementing  the  indices  associated 
with  the  respective  spatial  coordinate  directions  by  1,  for  example, 
u(i-l,j,k)  denotes  the  velocity  component  in  the  x-direction  associated 
with  the  left  cellular  face  of  the  (i, j ,k)~cell.  The  remainder  of  this 
section  is  devoted  to  presenting  each  of  the  finite  difference  forms  of 
the  governing  partial  differential  equations. 

III. 4.1  Conservation  of  Mass 

The  continuity  equation,  Eq  (13),  is  the  only  equation  to  retain 
its  conservative  form.  The  time  advanced  value  of  the  mixture  density 
for  any  cell  is  then  determined  from  the  following  difference  equation: 


Piti.k-Qi.i.k  +  (pu  >i!ik-  <pu 


At 


Ax 


+  ( pv  ^l.k  - ( pv  >£j-i,k 

Ay 


(  pv  )£},t  -  (  P» 


+ 


o 


(59) 


The  mixed  flux  expressions  (  Pu  ),  (  Pv  ),  and  (  Pw  )  are  differenced  in 
accordance  with  the  SOLA-ICE  algorithm.  As  an  example,  the  quantity 

vU+1 


(  Pu  ):£  j  ^  is  formulated  as  follows: 


<  pu  )^‘_k  .  +  )  "tj.k  + 


a  |  p 


n  /  ,,n+l  I  ,.n+l 

i+l.j,k  (  ui,j,k  I  ui,3fk 


I  )  + 


ci,J.k  (  u°tj,k  +  1  “”tj,k  1  > 


(60) 


The  remaining  flux  terms  of  Eq  (59)  are  differenced  similarly  and  index 
adjustments  follow  directly.  The  donor  cell  parameter,  a ,  governs  the 
amount  of  flux  advection  between  adjacent  cells;  it  is  further 
constrained  as  follows: 


(X  ®  \ 


0,  pure  spatially  centered  differences 


i ,  pure  donor  cell  differences 


The  donor  cell  terms  are  essentially  upwind  difference  foies.  The 
requirement  to  check  the  direction  of  the  flow  between  cells  is 
incorporated  in  the  program  without  adding  any  conditional  statements 
which  slow  the  cycle  computational  time.  The  method  as  utilized  in  the 
SOLA-ICE  algorithm  and  depicted  in  Eq  (60)  both  adds  and  subtracts  the 
absolute  value  of  the  time  advanced  velocity.  Dependent  upon  the 
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direction  of  the  fluid  flow,  the  appropriate  value  of  the  density  of  the 
fluid  is  provided  by  the  donor  cell. 


III. 4. 2  Conservation  of  Momentum 

The  three  momentum  conservation  equations  in  non-conservative, 

primitive  variable  form,  Eq  (23),  are  approximated  using  a  combination 

of  donor  cell  and  standard  difference  forms c  Again,  advection  terms  of 

the  form  U.  U.  .  and  U.  U.  .  use  the  modified  donor  cell  formulation, 
i  i.j  i 

whereas  the  remaining  terms  use  standard  upwind  or  central  difference 
forms  as  appropriate.  The  following  cellular  face  centered  flux 
expressions  of  the  velocity  gradients  are  provided  to  illustrate  the 
methodology  involved. 

If  the  x-direction  momentum  equation  is  arbitrarily  chosen,  the 


following  advection 

terms  occur: 

/u  5u\ 

n 

n 

»  /v  3u  \  ,  / 

n 

'w3u\ 

\  3x/ 

i.jfk 

V  3y/i,j,k  \ 

3z  /  i.j ,k 

These  terms , 

unlike 

those  encountered  in 

the  continuity  equation 

which 

was  in  conservative 

(divergence)  form,  are 

approximated  in  the  following 

manner : 

/u  3u\ 

\  3*/ 

n 

"  ui».1»k/ui+l,1,k  ”  u?-l,j,k\(  1  "  a  ) 

'  2Ax 

/ 

+  a  “  ui,i,k 

\  r ".j.t  - 1  ui,.i.k  1 1 
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(61) 

\  Ax 

/  \  2  / 

40 


I 


% 

1 

1 

ttH 

1 

Vi,j.k 

1 

g 

t 

& 

i.j.k 

1 

K 

1 

where  0  <  a  ^  1  . 

The  two  remaining  expressions, 

n  n 

/ v  JS  \  and  /  w  _9u  \ 

\  ty )  1,1, k  \  9z'i,j,k 

are  similarly  approximated,  however,  the  v  and  w  velocity  components  are 
not  defined  in  the  spatial  marching  direction  (in  this  case,  the  x- 
direction)  and  therefore  require  special  attention.  Only  the  u 
component  velocity  is  defined  in  the  x-direction  and  it  is  physically 
centered  on  both  the  left  and  right  faces  of  any  particular  cell. 
Consequently,  both  v  and  w  must  be  formed  from  the  averages  of  the  four 
v  and  w  neighboring  velocities  in  both  the  j -plane  and  k-plane 
respectively.  These  average  velocities,  v  and  w  are  defined  as  follows: 


<  ''i+i.i.k  +  ^l+i.i-i.k  *  vi.i.k  +  >  <62> 


4 


(63) 


With  these  definitions  of  vj  j  ^  and  w^  ^  ^  the  required  difference 
expressions  can  therefore  be  formed  analagously  to  Eq  (61),  that  is. 
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and, 
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(65) 


Eq  (65)  would  be  conceptually  valid  if  the  mesh  in  the  z-direction  was 
uniform,  however,  the  mesh  in  this  direction  varied  geometrically  as 
prescribed  by  Eq  (22).  The  proper  form  of  the  finite  difference 
approximation  for  the  advected  gradient  flow  in  the  z-direction  required 
replacement  of  the  uniform  mesh  spatially  centered  and  one-sided 
differences  by  the  non-uniform  mesh  formulations  of  section  III. 2  . 

Upon  defining  the  quantities  A z+  and  A z”  as  follows. 


zc(k+l)  -  zc(k) 
zcW  -  zr(k-l) 


(66) 


Eq  (65)  was  rewritten  in  the  following  form  which  is  reducible  to  that 
for  the  uniform  z-direction  mesh  when  K  is  identically  one. 
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where  Az+  and  Az~  are  defined  by  Eq  (66)  and  A  z^  is  the  vertical 
dimension  of  the  (i, j ,k)-cell.  The  variances  in  the  formulations  occur 
naturally  from  the  staggered  uniform-geometric  mesh  system  which  is 
utilized.  Care  must  be  exercised  in  forming  the  derivatives  in  the 
remaining  two  spatial  marching  directions.  Some  forms  are  identical  to 
those  listed  above,  whereas  others  are  different  and  require  additional 
inputs  into  the  numerical  simulation  model.  Incorporation  of  each  of 
the  above  finite  difference  approximations  and  the  remaining  required 
difference  forms  not  explicitly  presented  was  easily  accompished  for  the 
set  of  momentum  equations.  The  resultant  algebraic  expressions  were 
explicitly  solved  for  the  time  advanced  values  of  the  component  cell 
velocities  . 
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III. 4. 3  Conservation  of  Energy 

The  energy  equation,  Eq  (25),  was  approximated  by  its  finite 
difference  counterpart  in  the  same  manner  as  were  the  momentum 
equations.  In  this  case,  however,  the  specific  internal  energy  of  each 
cell  was  donor-cell  fluxed  as  were  the  advection  terms  of  the  momentum 
equations.  Terms  of  the  form, 


3l  +  u  3l  4-  v  3l  +  w  3l 
3t  3x  3y  3z 

were  approximated  as  follows  for  uniform  meshes  in  each  of  the  spatial 
directions  utilizing  the  shorthand  subscript  notation  whereby  omitted 
indices  are  held  constant: 


-  r°...k  +d-«)  |  ( uf 1  -k  u£j )  ( i;+1  - 1;_,  > 
At  2  2AX 


,n+l  .  „n+l  >  ,  Tn 


+  a 
2 


Tuf-1 

,  n+1  _ 

+  Ui-1 

1  u£+1 

+  »n+l 

r 

rTn  _  Tn 
ii+l  xi 

2 

2 

• 

Ax 

•  « 

u"+1  + 

«n+1  +  1 
ui-l  +  1 

u*+1  + 

un+l  | 1 
ui-l  1 

rTn  _  Tn 
xi  xi-l 

2 

2 

Ax 

+  (  1  -a  ) 


,,n+l  .  ,.n+l  Tn  _  Tn 

I vi  *  vi-i  \  (h+Llhzi.) 

\  2  /  l  24y  / 


+  a_ 
2 


ryf1  4-  -  |  v1 


,n+l  4.  ,.n+l 

i  vi-i 


In  -  Tn 
Ay 


r  vn+l  i  vn+l  .  |  „n+l  .  vn+l 

T.i  y.m  * 1  v.i  Y.i-i 


Tn  _  Tn 

1l._J.-i 

Ay 


45 


+  (  1  -  a  ) 


^  vk+1  +  vk-l  j  ^  Zk+1  ~  :k-l  j 


+  a 
2 


wk+1  +  «£li  -  i  wk+1  +  wk-l  i 


r  Tn  _  Tn 
Ik+1  xk 

A  z 


“f 1  +  "£-»  +  I  “f1  +  ”£{  I 


Tn  -  Tn 
\  ^k-l 

A  z 


(80) 


However,  since  a  geometrically  varying  mesh  was  utilized  in  the  z- 
direction,  the  spatially  centered  z-derivatives  must  be  approximated 
with  the  formulation  presented  in  section  7.11. 2,  and  the  increment  Az 
replaced  by  either  Az+  or  Az~  io  the  windward  differences  as 
appropriate. 

The  remaining  terms  of  the  energy  equation,  including  both 
components  c£  the  heat  flux,  were  approximated  in  the  same  manner  as 
were  the  diffusive  terms  of  the  momentum  equations.  The  conductive 
components  of  the  he.it  flux  terms  were  obtained  by  expanding  the 
tensorial  form  and  carrying  out  the  indicated  differentiations  and  then 
replacing  each  component  by  its  spatially  centered  finite  difference 
form.  Those  forms  in  the  z-direction  were  appropriately  modified  to 
account  for  the  geometric  mesh  variation.  The  diffusive  transport  of 
enthalpy  term, 

N 

(Jj^iff  k  *  hk  >’j  (81) 

k 


was  first  simplified  for  a  binary  gas  mixture  and  rewritten  as. 
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(  i  [  ^2  —  hj  ]  ),j  (3.2) 

where  ^  is  the  diffusive  flux  of  species  "1"  and  h^  and  are  the 

specific  enthalpies  of  species  "1"  and  species  ''2"  respectively.  The 
last  expression  was  then  expanded  utilizing  the  definition  of 
that  is, 


*dif  f  1  ”  "P®8”! 


3xj 

where  P  *  ^2  *  ^21  ^or  a  ^nary  8as  fixture  and  where  is  the  mass 
fraction  of  species  "1"  gas,  resulting  in  the  following  form  which  was 
then  approximated  by  finite  difference  expressions: 
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(83) 


The  derivatives  of  the  product  (pP  )  were  not  expanded  further,  but 
rather,  retained  in  this  form  since  the  gradients  of  this  quantity  would 
pose  less  difficulty  in  numerically  evaluating  them  as  would  evaluating 
each  of  the  expanded  gradient  terms  separately.  Additionally,  the 
spatially  centered  and  one-sided  derivative  terms  in  the  z-direction 
were  modified  accordingly  to  account  for  the  geometric  variation  of  the 
mesh  in  that  coordinate  direction. 


[yfiffizs- 
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III. A. 4  Conservation  of  Species 


The  conservation  of  species  partial  differential  equation  for  the 
ktn  species,  Eq  (22),  was  also  approximated  by  a  corresponding  finite 
difference  equation.  In  this  case,  the  differencing  was  formulated 
similarly  to  that  of  the  advection  terms  of  the  energy  equation 
utilizing  the  donor-cell  forms  developed  earlier  for  the  SOLA-ICE 
algorithm.  Rivard,  Butler,  and  Farmer  (Ref  46)  also  utilized  difference 
forms  for  the  species  conservation  equation  which  were  handled  similarly 
to  those  they  developed  for  their  energy  conservation  equation;  however, 
their  program  (RICE)  was  not  donor  cell  fluxed. 

^  The  general  form  of  this  conservation  equation  which  temporarily 
includes  the  source  terms  from  chemical  reactions  to  show  all 
contributory  effects  is  written  as  follows: 

(Pk),t  +  <pkUj)»j  “  mk,j  >’j  +  Sk  '-84) 


advection  diffusion  source 


The  actual  finite  difference  construction  to  represent  this  partial 
differential  equation  was  accomplished  in  two  steps  after  neglecting  any 
and  all  source  terms.  Firstly,  the  advective  fluxing  was  modeled  in  the 
same  manner  as  the  global  energy  advection  terms  employing  the  donor¬ 
cell  methodology.  Secondly,  the  diffusive  fluxes  were  modeled  similarly 
to  the  stress  tensor  terms  of  the  momentum  equations.  The  partial 
derivative  expansions  for  both  the  advective  and  diffusive  terms  are 
presented  respectively  as  follows: 

Advection  : 
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(85) 


(pk  uj  >m  “ 

P  k  f3u+3v+3w1+u3pk  +  v3pk+w3pk 
L  9  x  9y  3  z  j  3x  3y  3z 

Diffusion  : 


(pPmkij  ),j  - 

,3pS  .  +  ,&£.  /^\  +/3°I_  W3”k  \ 

'3  x  '  '  3  x  '  \3y/\3y/  '  3  z  M  3  z  ' 

+  pV  ["  3  +  3  2mk  +  3  2mk1  (86) 

l_3x2  3  y2  3  z 1  J 

As  in  the  global  energy  equation,  the  quantity  (p£>  )  was  treated  as  a 
property  of  the  fluid  and  difference  expressions  were  computed 
accordingly  for  both  the  uniform  and  geometrically  varying  meshes  in  the 
respective  spatial  coordinate  directions. 

III. 4. 5  Boundary  Conditions 

The  "no  slip",  continuative  outflow,  and  prescribed  inflow  boundary 
conditions  were  implemented  through  the  use  of  "dummy  cells"  as  opposed 
to  using  higher  order  difference  approximations.  The  "no-slip" 
condition  on  the  injection  surface  required  that  the  following  sub¬ 
conditions  be  strictly  enforced: 

(1) .  the  normal  velocity  vanish  on  the  surface,  and, 

(2) .  the  tangential  velocities  vanish  on  the  surface. 
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Implementation  of  the  boundary  conditions  associated  with  the  upper 
surface  of  the  computational  domain  was  accomplished  in  two  steps.  The 
first  step  covered  the  entire  upper  surface  with  a  "no-slip"  boundary 
which  was  physically  accomplished  by  the  following  mechanism: 

w(i,j ,kmax-l)  *0  ,  for  every  i,j 

u(i,j,kmax)  *  -u(i, j ,kmax-l)  ,  for  every  i,j  (87) 

v(i,j,kmax)  *  -v(i, j ,kmax-l)  .for  every  i.j 

Additionally,  the  variables,  P  and  I,  were  constrained  to  enforce  the 
adiabatic  wall  by  the  following  numerical  prescription: 

I(i»j,kmax)  *  I(i, j ,kmax-l)  ,  for  every  i.j 

p  (i,j,kmax)  =  P  (i,j,kmax-l)  ,  for  every  i.j  (88) 

The  second  step  redefined  the  variables  at  the  jet  orifice  to  conform  to 
the  prescribed  jet  inflow  values.  Each  of  the  affected  variables  was 
overwritten  with  the  requiste  exhaust  jet  data. 

The  windward  face  of  the  computational  domain  was  similarly 
manipulated  to  enforce  the  prescribed  boundary  defined  by  the  aircraft 
altitude  and  Mach  number.  The  parameters  in  the  (1,1 ,k)  and  (2,j,k) 

locations  were  enforced  to  equal  one  another  for  every  "j.k"  in  the 

domain  for  all  time.  This  numerical  prescription  embodies  the 
assumption  that  the  windward  face  is  sufficiently  far  upstream  to  be 
unaffected  by  the  flow  injection  and  interaction  downstream'-  In  all 


>'■*  <nff  v  Wp»w^&  ■  . 


cases  examined,  this  boundary  was  physically  placed  to  ensure  that  this 
assumption  was  valid.  The  ambient  conditions  enforced  on  this  Doundary 
were  obtained  from  the  atmosphere  subroutine  described  in  Appendix  A 
which  details  the  1962  Standard  atmosphere  properties  as  a  function  of 
aircraft  pressure  altitude. 

The  remaining  boundaries,  all  of  which  are  continuative  outflow, 
were  prescribed  initially  as  follows: 


Front  Lateral  Face  : 


v(i,l,k)  *  v(i,2,k)  ,  for  every  i,k 


u(i,l,k)  =  u(i,2,k)  ,  for  every  i,k 


w(i,l,k)  =  w(i,2,k)  ,  for  every  i,k 


P(i,l„k)  =  P(i,2,k)  ,  for  every  i,k 


where  P  is  any  property  of  the  fluid  mixture,  e.g.  I,  etc. 


Back  Lateral  Face  : 


v(i, jmax-l,k)  **  v(i, jmax-2,k)  ,  for  every  i,k 


u(i,jmax,k)  =  u(i, jmax-1 ,k)  ,  for  every  i,k 


w(i,jmax,k)  »  w(i, jmax-1, k)  ,  for  every  i,k 


P (i,jmaxsk)  «  P (i, jmax-1, k)  ,  for  every  i,k 


Lower  Face  : 
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w(i,j,l) 

u(i,j,l) 

v(i,j,l) 

P(i,j,D 

Leeward  Face  : 


-  w(i,j,2) 
“  u(i,j,2) 
»  v(i, j ,2) 

-  P(i,j,2) 


for  every  i,j 
for  every  i,j 
for  every  i»j 
for  every  i,j 


(91) 


u(imax-l,j ,k)  =  u(imax-2, j ,k)  ,  for  every  j,k 

v(imax,j,k)  ■  v(imax-l, j ,k)  ,  for  every  j,k 

w(imax,j,k)  *  w(imax-l, j ,k)  ,  for  every  j,k 

P(imax,j,k)  *  P  (imax-*l,j,k)  ,  for  every  j,k  (92) 

After  the  solution  advanced  through  both  the  pressure  iteration  phase  of 
each  computational  time  cycle  and  the  subsequent  solution  of  the 
momentum  equations,  the  continuative  boundary  conditions  were  imposed. 

III. 5  Stability  Requirements 

To  ensure  the  stability  of  the  numerical  algorithm,  the  three 
following  constraints  were  strictly  enforced  during  each  time  step  and 
for  every  cell  within  the  computational  domain: 

(1) .  The  fluid  must  not  travel  a  distance 

greater  than  any  cell  dimension, 

(2) .  The  fluid  momentum  must  not  diffuse  a  distance 

greater  than  any  cell  dimension,  and, 
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(3).  The  fluid  energy  must  aot  diffuse  a  distance 
greater  than  any  cell  d-'.’msion. 


Explicit  mathematical  ex  vessions  for  these  limits,  with  the  exception 
of  the  mass  transport  limit,  are  not  possible  to  formulate  for  the 
coupled  non-linear  differential  equations  of  the  model.  However,  when 
the  momentum  and  energy  equations  were  approximated  by  the  corresponding 
pure  diffusion  counterparts,  Cloutman  et  al  (Ref  44)  provided  two- 
dimensional  mean  flow  expressions  derived  from  von  Neumann  stability 
analyses  which  were  easily  extended  to  fully  three-dimensional  turbulent 
flows  comprised  of  one  fluid  species.  These  expressions  which  are  meant 
to  serve  only  as  guides  and  not  exact  mathematical  prescriptions  are 
presented  as  follows: 

x-momentum  : 


f- 

'  1  u  | 

+  |v| 

+  |w| 

[2 

.Ax 

max  A  y 

max  A  z 

max  J 

P  eff  T  4  1  +  1  +  1  1  1 At  <  1  (93) 

P  L  3  Ax7  Ay2  aT2  J  j  1 


/-momentum 


a  J*  I u 1 1  +  |v|  +  |w|  I  "I  + 

2  L  Ax  I  max  Ay  max  Az  |maxj 

U  eff  f  1  +  4  1  +  l  ]  ]  At  <  1  (94) 

T - “LZ?  IA~7  Al2]/  2 
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z-momentum  : 


fa  r  |u| 

+ 

|v| 

+  |w| 

1 

{2  L  zr 

max 

Ay  max 

A  z 

max  J 

y  eff 
P 
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+  1  4 

Ay2  3iz2_ 

]  At 
} 

energy  : 


(95) 
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(96) 


As  the  direct  result  of  numerical  experimentation  undertaken  to  verify 
the  utility  of  the  above  expressions,-  it  was  observed  that  a  time 
increment  less  than  or  equal  to  one-fourth  the  value  obtained  by 
satisfying  the  cellular  mass  transport  constraint,  that  is. 


t  <  min  J^x  ,  Ay  ,  A  z  1  (97) 
t  M  |v|  |w|  J 

provided  an  effective  means  to  determine  the  proper  temporal  increment 
to  ensure  numerical  stability.  If  this  time  step  when  substituted  in 
the  stability  analysis  equations  defined  above  results  in  a  donor-cell 
parameter  value  which  exceeds  its  constraint,  that  is  &  >  1,  then  the 
temporal  increment  must  be  decreased  such  that  the  donor-cell,  momentum, 
and  energy  constraints  are  satisfied  simultaneously. 
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Because  of  the  violent  nature  of  the  high-speed  injection  and 
associated  viscous  flow  interactions,  the  coarseness  of  the  vertical 
grid  directly  affected  the  stability  of  the  numerical  solution.  For 
high  Mach  number  crossflows  (0.7)  with  respect  to  the  jet  injection  Mach 
number  (0.24),  a  cell  stretch  factor  up  to  and  through  1.2  provided 
accurate  and  stable  solutions.  However,  for  low  Mach  number  crossflows 
(0.125)  with  respect  to  the  jet  injection  Mach  number  (0.5),  a  stretch 
factor  of  1.1  was  too  large.  In  the  latter  case,  a  stretch  factor  of 
1.05  proved  to  be  small  enough  to  capture  the  essential  near  field 
mixing  and  to  guarantee  both  the  accuracy  and  stability  of  the  numerical 
solution.  From  this  discussion,  it  should  be  obvious  that  each 
injection  problem  must  be  carefully  analyzed  on  a  case  by  case  basis  and 
that  those  parameters  which  can  be  altered  to  affect  stability  must  be 
"tuned"  for  each  problem. 

III. 6  Convergence  Criteria 

The  cost  of  solving  three-dimensional  fluid  flow  problems  is 
exhorbitant  and  consequently  suitable  criteria  are  required  to  determine 
when  to  cease  numerical  computations.  Simply  restated,  there  exists  a 
point  in  time  beyond  which  the  marginal  returns  from  additional 
computations  do  not  offset  the  increased  cost  of  acquiring  additional 
information.  Many  authors  have  used  the  criterion  of  the  asymptotic 
approach  to  "steady  state"  of  the  fluid  density  from  one  time  step  to 
the  next  whereby  the  difference  in  density  is  less  than  some  arbitrarily 
prescribed  increment.  However,  the  jet-crossflow  interaction  problem 
for  real  gases  required  that  some  criterion  be  established  to  determine 
not  only  when  the  change  in  density  was  less  than  some  prescribed  value, 
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but  also,  when  the  flow  field  was  firmly  established  to  accurately 
permit  computation  of  the  jet  trajectory.  The  criterion  chosen  to 
accomplish  both  of  these  requirements  was  the  asymptotic  change  in 
kinetic  energy  of  the  computational  domain  from  one  time  step  to  the 
next  being  less  than  some  prescribed  value.  This  criterion  is 
mathematically  expressed  as 


n*At  (n+1)  *A  t 

(  Kinetic  Energy  )  -  (  Kinetic  Energy  ) 


(  Kinetic  Energy  ) 


where  e  is  an  arbitrarily  small  positive  value. 

This  criterion  alone  is  insufficient  to  determine  the  "steady- 
state"  solution  to  the  problem  because  of  both  the  oscillatory  behavior 
of  the  kinetic  energy  and  the  small  temporal  increments  utilized.  Local 
maxima  and  minima  erroneously  would  indicate  convergence  by  satisfying 
Eq  (98)  when  obviously  it  did  not  exist.  If  termination  of  the 
computations  occurred  too  soon  based  upon  the  information  contained 
within  the  numerical  transient  solutions,  meaningful  information 
concerning  the  flowfield  and  its  associated  properties  would  not  be 
available.  However,  by  specifying  some  meaningful  minimum  time  to 
terminate  computations  based  upon  the  dynamics  of  the  simulation  model, 
this  pitfall  was  easily  avoided. 

The  time  required  for  a  disturbance  introduced  into  the 
computational  domain  to  decay  to  negligible  levels  is  directly  related 
to  the  characteristic  time  of  the  computational  domain  which  is  defined 
as  follows; 
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0. 

where  L  is  the  length  of  the  computational  domain,  and, 
where  is  the  magnitude  of  the  freestream  velocity. 


(99) 


It  was  observed  that  depending  upon  the  coarseness  of  the  computational 
grid,  the  measure  of  convergence  specified  by  Eq  (98)  was  obtained  by  at 
least  one  characteristic  time.  For  the  coarse  grids,  the  required 
convergence  level  was  attained  within  two  characteristic  times. 
Consequently,  when  Eq  (93)  was  utilized  in  conjunction  with  some 
multiple  of  t^  as  the  cutoff  time,  termination  of  the  computations 
was  effectively  accomplished  with  essentially  no  additional  information 
to  be  gained  from  additional  computational  effort.  These  criteria  were 
employed  on  numerical  experiments  involving  variable  property  air-air, 
helium-air,  and  gas  mixture-air  jet  injection  problems  at  elevated 
temperatures  (1500°K)  with  outstanding  success. 

III. 7  Summary 

The  computational  domain  defined  in  Section  II  was  divided  into  a 
large,  but  finite,  number  of  rectangular  cells,  the  vertical  dimensions 
of  which  ware  prescribed  to  vary  geometrically  from  the  injection 
surface.  The  unique  finite  difference  forms  required  to  approximate  the 
terms  of  the  non-linear  governing  partial  differential  equations  as 
implemented  by  the  S0La~ICE  algorithm  (Ref  44)  and  as  modified  to 
account  for  the  non-uniform  mesh  in  the  z-dfrection  were  presented.  The 
required  equivalent  spatially  "centered”  and  one-sided  difference  forms 
for  both  the  first  and  second  partial  derivatives  were  derived  by 


utilizing  Taylor  series  expansions  of  several  variables.  The  staggered 
mesh  employed  by  the  algorithm  presented  additional  complexity  to  the 
modeling  as  the  spatial  locations  of  variables  with  the  same  indices 
were  not  coincident;  consequently,  the  necessary  modifications  to  the 
finite  difference  approximations  utilized  by  the  algorithm  were 
presented.  The  boundary  conditions  were  easily  implemented  through  the 
use  of  "dummy"  cells  which  surrounded  the  computational  domain  as 
opposed  to  the  use  of  higher  order,  multi-point  difference 
approximations  thereby  increasing  the  computational  speed  by  the 
reduction  in  complexity  afforded  by  this  method. 

Since  a  modified  donor-cell  formulation  was  used,  the  algorithm 
required  the  strict  enforcement  of  the  constraints  that  no  cellular 
mass,  momentum,  or  energy  diffuse  a  distance  greater  than  any  cell 
dimension  during  any  given  time  step  for  each  and  every  cell  within  the 
computational  domain-.  Explicit  numerical  stability  expressions  derived 
from  von  Neumann  stability  analyses  of  the  two-dimensional  pure 
diffusion  counterparts  of  the  momentum  and  energy  equations  (Ref  44,  Eqs 
(52)-(54)  )  were  expanded  for  three-dimensional  turbulent  flows.  The 
rationale  for  choosing  the  extremely  successful  convergence  criterion 
for  this  problem,  the  asymptotic  approach  to  "steady-state"  of  the 
change  in  kinetic  energy  of  the  computational  domain  from  one  time  step 
to  the  next,  in  conjunction  with  a  computational  cutoff  time  based  on 
some  multiple  of  the  domain  characteristic  time  was  also  detailed.  The 
results  of  the  research  of  the  jet-crossflow  interaction  problem  for 
real  gases  at  elevated  temperatures  are  presented  in  the  next  section. 
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IV.  Results  of  the  Numerical  Simulation 

This  section  describes  the  numerical  solution  of  the  jet  injection- 
freestream  interaction  for  the  chemical  laser  exhaust  gas  mixture 
described  in  Table  I  as  well  as  for  both  pure  helium  and  air  jets 
constrained  by  the  same  parameters  defining  the  problem,  that  is,  the 
Mach  numbers  of  both  the  jet  and  the  freestream,  the  ratio  of  jet  to 
freestream  dynamic  pressure,  and  the  jet  injection  temperature  were 
fixed.  Numerical  analyses  demonstrated  that  the  trajectories  and 
general  characteristics  of  the  flowfields  of  the  jet-crossflow 
interaction  were  essentially  unaffected  by  the  turbulent  dif fusivities 
associated  with  the  thermal  conductivity  and  the  mass  transfer  binary 
diffusion  coefficient.  Consequently,  the  isopleths  from  the  simulations 
with  molecular  values  for  thermal  conductivity  and  binnrv  diffusion 
coefficient  for  the  AFWL  gas  mixture,  helium  and  air  injection  studies 
were  presented  as  the  trajectories  were  more  easily  perceived  as  well  as 
numerically  obtained.  Section  IV. 5  contains  the  discussion  of  the 
|  effects  induced  on  the  heat  transfer  to  the  injection  surface  by  the 

incorporation  of  the  energy  transport  diffusivities  associated  with  the 
turbulent  nature  of  the  flow  interaction.  To  verify  the  methodology  of 
the  algorithm,  numerical  simulations  were  conducted  for  an  air  jet 
injected  vertically  upwards  into  a  subsonic  crossflow  from  a  rectangular 
jet  orifice  with  an  aspect  ratio  of  4.0  and  with  a  jet  to  freestream 
velocity  (Mach)  ratio  of  4.0  (equivalently,  a  Q  ratio  of  16)  for  which 
experimental  flowfield  data  existed  and  from  which  direct  comparisons 
could  then  be  made. 

IV. 1  Verification  of  the  Numerical  Algorithm 


'tsx 
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The  experimental  normally  upwards  injection  tests  conducted  by 
Weston  and  Thames  (Ref  26)  at  NASA/Langley  for  an  aspect  ratio  4.0  jet 
with  the  major  axis  of  the  orifice  aligned  with  the  streamwise  flow  was 
numerically  simulated  as  indicated  in  Figure  6.  The  boundary  conditions 
for  the  computational  domain  were  identical  to  those  already  discussed 
for  the  downwardly  ejecting  jet  except  for  the  fact  that  the  upper  and 
lower  surfaces  were  interchanged  such  that  the  adiabatic  "no-slip" 
injection  surface  was  located  on  the  bottom  surface  through  which  the 
jet  was  forcibly  constrained  to  be  ejected  into  the  Mach  0.125  crossflow 
at  essentially  freestream  (tunnel)  static  temperature  and  pressure.  The 
particular  test  runs  chosen  were  for  a  jet  to  freestream  velocity  ratio, 
R,  equal  to  4.0  and  for  a  jet  injection  angle  of  90  degrees.  These  data 
were  for  the  lowest  values  of  R  experimentally  investigated  and  were 
therefore  more  nearly  representative  of  the  real  gas  injection  problem 
which  was  to  be  simulated. 

The  actual  tests  were  accomplished  over  a  time  span  of  several 
months,  consequently,  test  conditions  varied  with  the  climatic 
conditions  prevalent  during  each  test.  The  numerical  simulation 
mathematically  replicated  the  injection  tests  from  the  reflection  plane 
injection  surface  of  dimension  12De  by  24Dg  where  Dg  is  the  equivalent 
diameter  of  the  rectangular  jet.  Additionally,  average  values  of  the 
tunnel  static  temperature  (281.76  °K)  and  pressure  (2130  lbg/ft^)  were 
utilized  while  constraining  the  crossflow  (tunnel)  velocity  to  be  138 
feet  per  second.  Actual  experimental  flowfield  data  were  provided  by 
Weston  to  facilitate  direct  numerical  comparison  with  the  mathematical 
simulation. 


IV . 1 . 1  Numerical  Solution  :  Convergence 


The  experimental  tests  for  R=4  were  numerically  simulated  using 
both  a  coarse  mesh  of  29  x  15  x  17  cells  in  the  x,  y,  and  z  spatial 
directions  respectively  on  the  CDC  Cyber  175  at  Wright-Patterson  AFB  and 
a  more  refined  mesh  of  29  x  29  x  26  cells  in  the  x,  y,  and  z  directions 
respectively  on  the  CDC  7600  at  the  Air  Force  Weapons  Laboratory, 
Kirtland  AFB,  New  Mexico.  Because  the  experimental  jet  was  both  narrow 
and  small,  the  jet  orifice  was  modeled  by  2  x  1  cells  in  the  x  and  y 
spatial  directions  respectively  in  the  coarse  and  refined  mesh 
simulations.  Additionally,  in  each  case,  a  geometric  stretch  factor,  K, 
of  1.05  was  employed  to  expand  the  vertical  cell  dimension  from  the 
injection  surface;  however,  the  initial  cell  vertical  dimension  for  the 
refined  mesh  was  one-half  that  of  the  corresponding  value  of  the  coarse 
mesh.  Only  the  fine  mesh  simulation  replicated  the  actual  experimental 
injection  plane  dimensions. 

As  alluded  to  in  Section  III. 5,  the  numerical  stability  was 
directly  affected  by  the  magnitude  of  the  geometric  stretch  factor.  For 
the  coarse  grid  check  case,  a  K  value  of  1.1  resulted  in  a  solution 
which  initially  appeared  to  be  periodic  with  time.  The  convergence  plot 
of  this  case  is  presented  in  Figure  7.  After  approximately  five  (5) 
characteristic  times,  the  periodicity  of  the  solution  was  firmly 
established.  However  in  the  process  of  verifying  and  analyzing  the 
numerical  solution,  another  coarse  grid  test  simulation  was  made  with  a 
K  value  equal  to  1.05.  As  is  evidenced  by  Figure  8,  the  "periodic" 
solution  which  was  apparent  for  the  larger  K  factor  did  not  recur,  but 
rather,  the  asymptotic  approach  to  steady  state  resulted  after 
approximately  two  (2)  characteristic  times.  When  the  mesh  was  further 
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refined  by  both  halving  the  initial  cell  vertical  dimension  and 
increasing  the  number  of  cells  within  the  computational  domain,  the 
solution  exhibited  an  acceptable  level  of  convergence  within 
approximately  one  (1)  characteristic  time. 

IV . 1 . 2  Numerical  Solution  :  Results 

The  numerical  simulation  of  Weston's  experimental  test  for  a 
normally  injected  air  jet  with  R=4  produced  results  which  were  in 
agreement  with  those  observed  in  the  wind  tunnel.  As  the  jet  penetrated 
the  crossflow,  it  was  deflected  downstream  and  in  the  process,  a  pair  of 
counter-rotating  vortices  were  formed  which  also  were  swept  downstream. 
Figures  10a  through  lOe  depict  the  numerically  computed  flowfield  in 
adjacent  y-planar  cuts  through  the  computational  domain  where  Figure  10a 
is  situated  two  cells  to  the  left  of  the  jet  symmetry  plane.  As  the 
sweep  is  made  from  this  plane,  through  the  symmetry  plane,  and  finally 
concluding  with  the  plane  two  cells  to  the  right  of  the  symmetry  plane 
(Figure  lOe),  the  physical  nature  of  the  jet-crossflow  interaction  is 
readily  observed.  Behind  the  jet,  a  wake  region  with  a  small 
recirculation  zone  is  observed  which  entrains  low  velocity  air  near  the 
reflection  plane  into  the  jet  plume.  As  the  lateral  distance  from  the 
sides  of  the  jet  is  increased,  the  flowfield  demonstrates  less  influence 
from  the  relatively  high  speed  jet  injection  and  the  velocity  tends  to 
return  to  frees tream  values. 

As  Weston's  data  are  the  only  known  experimental  flowfield  data  for 
rectangular  jets  being  injected  into  crossflows,  it  was  imperative  that 
specific  data  be  compared  to  ascertain  the  accuracy  of  the  numerical 
solution  of  the  jet-crossflow  interaction.  The  symmetry  plane  was 
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chosen  because  of  the  additional  analysis  provided  by  Weston  and  Thames 
in  deriving  empirical  least  squares  determined  expressions  for  the  jet 
trajectory  (based  on  maximum  total  pressure)  as  well  as  for  the  locus  of 
centers  for  the  counter-rotating  vortex  pair  (Ref  26,  Eq  1  and  Table  2). 
Figure  11  is  a  vector  plot  of  experimental  wind  tunnel  data  (solid 
vectors)  and  flowfield  data  from  the  refined  mesh  numerical  simulation 
(dotted  vectors)  for  the  jet  injection  symmetry  plane.  Since  the 
spatial  locations  of  the  test  data  were  not  coincidental  with  those  of 
the  simulation,  a  two-dimensional  linear  interpolation  routine  was  used 
to  acquire  simulation  data  at  the  requisite  experimental  test  locations. 
Excellent  agreement  was  obtained  for  data  along  the  jet  trajectory  both 
ir>  angular  orientation  and  in  magnitude.  It  was  disconcerting  to 
observe  that  the  empirical  velocity  centerline  curve  also  plotted  in 
Figure  II  demonstrated  approximately  the  same  slope  as  the  velocity 
field  generated  by  the  numerical  solution  while  the  experimental  data 
upon  which  the  empirical  curve  was  based  showed  marked  deviation. 
However,  the  least  squares  analysis  performed  by  Weston  and  Thames  in 
deriving  the  trajectory  equation  was  based  upon  the  flowfield  data  of 
the  exhaust  plume  in  the  symmetry  plane  as  well  as  adjacent  lateral 
planes.  Consequently,  the  empirical  expression  more  closely  portrays 
the  actual  macroscopic  plume  behavior.  The  experimental  values  of  the 
velocities  obtained  by  the  first  probe  on  the  rake  (highest  Z/Dg  value) 
had  the  greatest  distortion  induced  measurement  error.  Other  near  field 
data  did  not  offer  the  same  degree  of  agreement  as  did  either  the  far 
field  or  jet  trajectory  data.  The  empirically  derived  curve  for  the 
locus  of  centers  of  the  vortex  pair  is  also  presented  within  the  same 


figure.  It  Is  readily  apparent  that  the  data  acquired  during  the 
experimental  phase  were  more  heavily  weighted  in  favor  of  determining 
the  structure  of  the  vortex  pair  than  determining  the  jet  centerline 
based  on  maximum  velocity.  Because  of  the  paucity  of  experimental 
flowfield  data  on  or  near  the  jet  trajectory  in  conjunction  with  the 
probe  errors  induced  by  the  large  gradients  over  the  probe  orifice  in 
the  near  field  region,  the  experimental  data  immediately  above  and  below 
the  empirical  curve  in  Figure  11  were  quantitatively  compared  both  in 
magnitude  and  flow  orientation  with  the  corresponding  data  from  the 
simulation.  The  magnitudes  of  the  velocities  in  the  streamwise 
direction  above  and  below  the  trajectory  in  Region  I  of  the  figure  were 
underpredicted  by  the  simulation  by  7.4  and  12  per  cent  respectively. 
The  local  flow  inclination  angles  from  the  simulation  were  7.9  and  9.5 
degrees  respectively;  the  corresponding  experimental  angles  were  -0.8 
and  1.3  degrees.  The  inclination  of  the  empirical  curve  was  19.3 
degrees  at  the  intersection  of  the  rake  probe  axis  with  the  trajectory. 
Similar  data  for  Region  II  of  the  sa~v.  figure  shows  that  the  simulation 
underpredicted  the  magnitudes  of  the  streamwise  velocities  by  2.3  and 
8.4  per  cent  for  the  upper  and  lower  data  respectively.  The 
corresponding  local  angular  orientations  varied  between  5.9  and  7.3 
degrees  from  the  simulation  versus  the  experimentally  obtained  -2.4  and 
1.0  degrees  while  the  inclination  of  the  empirical  curve  is 
approximately  16.9  degrees  at  the  intersection  of  the  trajectory  and  the 
rake  axis.  In  both  regions  along  the  empirically  derived  curve  for  the 
jet  trajectory,  it  is  readily  observed  that  the  numerical  simulation 
more  closely  replicates  the  trends  exhibited  in  the  curve  than  do  the 
experimental  data  for  the  symmetry  plane. 


To  demonstrate  the  effects  of  the  mesh  cell  sizing  on  the  numerical 
resolution.  Figure  12  is  presented  which  depicts  the  flowfield  in  the 
symmetry  plane  from  the  coarse  grid  simulation.  The  same  general  flow 
pattern  as  in  Figure  10c  was  discerned  but  not  to  the  same  level  of 
accuracy.  However,  the  same  small  recirculation  zone  entrainment  of  low 
velocity  air  near  the  injection  plane  was  still  observed  as  well  as  was 
the  same  wake  behavior  of  the  jet-crossflow  interaction  downstream  of 
the  jet. 

Weston  and  Thames  presented  pressuj.^  coefficient  (Cp)  data  for  the 
rectangular  jet  aligned  with  the  flow  (Ref  26,  Figure  6b).  However, 
their  presentation  was  for  the  ACp  derived  by  subtracting  Cp  data  with 
the  jet  off  from  the  corresponding  values  with  the  jet  on.  Despite  the 
fact  that  the  numerical  simulation  was  not  optimized  to  obtain  surface 
pressure  distributions,  that  is,  the  grid  was  still  too  coarse  even  in 
the  refined  mesh  case,  results  were  obtained  which  were  nearly  in 
agreement  with  experimental  surface  pressure  measurements.  Figure  13a 
depicts  Cp  contour  data  derived  from  the  mathematical  simulation.  A 
large  wake  region  (negative  Cp)  exists  laterally  and  downstream  of  the 
jet  while  a  small  positive  pressure  region  exists  immediately  upstream 
of  the  jet.  The  vortex  structure  is  also  observed  from  this  figure  as 
well  as  is  the  characteristic  kidney  shape  of  each  Cp  contour. 

Figure  13b  shows  the  CD  contour  data  for  the  wake  region  which 
Weston  and  Thames  included  in  Reference  26.  The  Cp  values  which 
they  subtracted  out  in  their  presentation  varied  between  +0.03  with  two 
anomalous  variations  of  magnitude  0.05  in  Cp.  Essentially,  the  data 
contained  in  their  presentation  are  equivalent  to  the  data  in  Figure 
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13a;  however,  the  ambient  temperature  and  pressure  at  which  the 
experimental  data  were  obtained  are  not  known.  Both  figures  demonstrate 
similar  behavior  with  the  experimental  values  exhibiting  more  pronounced 
and  intense  variations.  To  quantify  the  comparison  between  the 
experimental  and  analytical  results,  Cp  data  for  a  streamwise  station 
X/Dg  =  0.0  which  were  nrovided  by  Weston  were  depicted  in  Figure  13c 
with  the  corresponding  data  derived  from  the  simulation.  This  figure 
shows  what  appears  to  be  a  constant  difference  acting  in  conjunction 
with  a  rotation  in  the  data.  Because  of  the  velocity  averaging  scheme 
employed  in  the  simulation,  such  an  aberration  in  the  data  is  entirely 
consistent  with  the  coarsness  of  the  mesh  which  was  utilized.  As  the 
lateral  distance  increases,  both  sets  of  data  demonstrate  the  same 
trend. 

Another  similar  experiment  conducted  by  D.  J.  Peake  (Ref  47) 
experimentally  measured  the  surface  pressure  distribution  on  the  bottom 
surface  of  a  wing  in  which  an  aspect  ratio  4.0  rectangular  jet  was 
embedded  with  the  major  axis  aligned  with  the  streamwise  flow.  The  R 
values  for  his  tests  were  fortuitously  4.0  and  4.1;  the  crossflow 
velocities  were  60  and  120  feet  per  second,  respectively.  Figures  14a 
and  14b  were  reproduced  from  Peake's  work  (Ref  47,  Figures  10a  and  lOe) 
and  depict  the  pressure  coefficient  contours  experimentally  determined 
on  the  flat  region  of  the  lower  wing  surface.  The  quantity  which  Peake 
used  to  nondimensionalize  the  spatial  locations  was  not  the  equivalent 
diameter  of  the  rectangular  orifice  but  rather  was  the  square  root  of 
the  orifice  area. 

The  flow  over  the  wing  surface  in  which  Peake  embedded  the 
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rectangular  jet  was  affected  by  the  upstream  planform  geometry,  the 
effects  of  which  became  more  pronounced  as  the  crossflow  velocity  was 
increased  while  maintaining  the  R  value.  Because  the  crossflowing 
airstream  was  subsonic  (60  and  120  feet  per  second),  distortions  of  the 
pressure  coefficient  contours  induced  by  the  upstream  and  finite  edge 
conditions  associated  with  the  experimaental  apparatus  negated 
meaningful  comparisons  between  the  experimental  and  analytical  data  for 
regions  outside  the  near  wake  zone.  To  that  end,  comparisons  were  made 
for  X/De  and  Z/Dg  pairs  (2.25,1.1)  and  (2.75,1.0)  for  which  the 
analytically  determined  values  were  0.25  and  0.20  respectively  (see 
Figure  13).  At  the  corresponding  locations  in  Peake's  non-dimer.sional 
terminology,  (2.54,1.24)  and  (3.10,1.13)  respectively,  it  was  observed 
from  Figure  14a  that  for  a  crossflow  velocity  of  60  feet  per  second  the 
corresponding  values  of  the  pressure  coefficient  were  0.25  (nominal 
difference)  and  0.20  (nominal  difference).  Whereas  in  Figure  14b  at  the 
same  prescribed  coordinates,  the  corresponding  values  of  Cp  were  0.30 
and  0.24  with  associated  differences  of  +20  per  cent  each.  The  larger 
discrepancies  accompanying  the  increase  in  crossflow  velocity  (idO  feet 
per  second)  can  not  be  attributed  only  to  the  slight  increase  in  R  (4.1 
vice  4.0)  although  there  is  necessarily  some  tangible  effect.  The 
effect  of  the  flat  plate  injection  surface  on  the  flowfield  as  pointed 
out  by  Weston  and  Thames  (Ref  26)  is  not  as  benign  as  previously 
thought.  Therefore,  a  modified  wing  surface  with  only  a  small  region  of 
negligible  curvature  is  even  more  insidious.  Additionally,  Fearn  and 
Weston  (Ref  48)  have  shown  that  errors  in  comparisons  between 
corresponding  Cp  contours  on  the  surface  of  a  flat  plate  in  regions  of 
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low  pressure  gradients  for  two  identical  experiments  at  the  same  set  of 
conditions  can  very  easily  vary  in  spatial  orientation  by  as  much  as  an 
equivalent  diameter.  Consequently,  the  overall  agreement  oetween  the 
simulation  and  the  experiment  is  even  further  enhanced. 

Perhaps  the  most  information  regarding  the  surface  pressure 
distribution  was  obtained  from  Figure  15  which  presents  a  three- 
dimensional  0^  surface  plot  for  the  injection  plane  and  the  associated 
c-julour  plot  of  it  p-_ ojected  onto  the  x-y  plane  beneath  it.  The  surface 
plot  demonstrates  the  typical  flat  plat  pressure  distribution  which  is 
characterized  by  a  pressure  rise  at  the  leading  edge  which  decays 
rapidly  to  a  null  value.  Just  ahead  of  the  jet,  a  positive  pressure 
region  exists  which  abruptly  becomes  negative  aft  of  the  jet  orifice.  A 
pressure  disturbance  issues  forth  laterally  from  the  sides  of  the  jet 
which  rapidly  decays  prior  to  reaching  the  lateral  plate  boundaries. 
Also  observable  from  the  figure  is  a  positive  pressure  region  aft  of  the 
jet  which  forms  a  "spine-like"  surface  which  gradually  dissipates  over 
the  entire  'eraaining  length  of  the  injection  plane.  This  turbulent 
region  was  also  observed  by  Weston  in  the  wind  tunnel  (Ref  49)  but 
unfortunately  «.r  was  measured  qualitatively.  In  the  contour 

plane,  a  contour  of  this  positive  Cp  zone  is  displayed  as  {••re  the 
characteristic  kidney  shapes  of  the  wake  region  and  the  positive  C 
region  contours  ahead  of  the  jet. 

IV. 2  Airborne  Chemical  Laser  Exhaust 

The  analysis  of  the  high  temperature  (1500  °K)  exhaust  gases 
generated  by  an  airborne  chemical  laser  and  subsequently  discharged  from 
an  aspect  ratio  1.75  rectangular  jet  orifice  embedded  in  the  aft  lower 
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surface  of  a  wide-bodied  air  vehicle  fuselage  into  the  high  subsonic 
(M=0.7)  ambient  crossflow  at  a  pressure  altitude  of  11  kilometers 
(tropopausa)  is  presented  within  this  subsection.  The  computational 
domain  and  the  requisite  finite  difference  formulations  of  both  the 
governing  equations  »nd  the  boundary  conditions  have  been  described  in 
Sections  II  and  III.  The  inherent  differences  between  this  study  and 
that  previously  conducted  to  verify  the  algorithm,  that  is,  the  vertical 
isothermal  air  injection  problem  for  R=4,  are  significant.  This  case 
(and  the  succeeding  ones)  involved  a  high  temperature  subsonic 
rectangular  jet  being  injected  into  a  cold  crossflowing  gas  not 
necessarily  comprised  of  the  same  species  as  that  of  the  injectant.  The 
Q  ratio  was  also  markedly  different;  Q  for  this  and  the  subsequent 
studies  was  constrained  to  the  value  of  0.15  while  the  corresponding  Q 
value  for  the  Weston  experiment  was  16.  Additionally,  the  scales  of 
both  problems  differed  appreciably.  For  example,  the  boundary  layer 
thickness  upstream  of  the  high  temperature  jet  was  approximately  equal 
to  D/14  where  D  is  the  streamwise  dimension  of  the  rectangular  jet;  in 
Weston's  experi.raent,  the  corresponding  upstream  boundary  layer  thickness 
was  approximately  D/36  where  D  again  is  the  streamwise  dimension  of  the 
rectangular  jet.  However,  the  boundary  layer  thickness  associated  with 
the  simulations  for  hot  gas  -  cri  injectants  was  two  orders  of  magnitude 
greater  than  that  of  the  isothermal  case.  In  dimensional  terms,  the 
actual  boundary  layer  thickress  was  approximately  0.3  meters  for  the 
airborne  chemical  laser  simulations. 

IV. 2 - l  Numerical  Solution  :  Convergence 

The  numerical  simulation  of  the  exhaust  jet-crossflow  interaction 
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was  accomplished  for  both  coarse  and  refined  subdivisions  of  the 
computational  domain.  The  coarse  mesh  simulation  conducted  on  the  CDC 
Cyber  175  at  Wright-Patterson  AFB  utilized  29  x  13  x  13  cells  in  the  x» 
y,  and  z  spatial  coordinate  directions  while  the  refined  mesh  simulation 
conducted  on  the  CDC  7600  at  Kirtland  AFB  used  42  x  21  x  23  cells  in  the 
corresponding  spatial  directions.  The  jet  orifice  was  modeled  by  4  x  2 
and  5x3  cells  in  the  x  and  y  directions  for  the  coarse  and  refined 
meshes  respectively.  The  K  factors  controlling  the  geometrically 
varying  z-dimcnsion  of  each  cell  from  the  injection  surface  were  1.20 
and  1.09  for  the  coarse  and  refined  mesh  simulations.  Despite  the  large 
value  of  K  utilized  by  the  coarse  mesh  simulation  which  was  required  to 
ensure  that  the  near  field  mixing  processes  were  essentially  captured 
within  the  computational  domain,  the  criterion  to  establish  convergence 
(Figure  16)  exhibited  both  numerical  stability  and  the  associated 
approach  to  steady-state  after  approximately  two  (2)  characteristic 
times.  The  corresponding  data  for  the  refined  mesh  (Figure  17) 
exhibited  an  acceptable  level  of  convergence  within  one  (1) 
characteristic  time.  As  is  evidenced  by  Figure  17,  the  convergence  of 
the  solution  was  accelerated  by  the  modification  of  the  three  factors 
which  governed  both  the  stability  and  accuracy  of  the  algorithm,  and 
consequently,  of  the  numerical  solution.  From  time  t=*0  to  approximately 
t=0.03  seconds,  the  donor-cell  parameter  was  equal  to  C. 7,  the  over- 
relaxation  factor  was  1.4,  and  the  temporal  increment  was  .-.qual  to 
0.00015  seconds.  Because  the  magnitude  of  the  change  in  domain  kinetic 
energy  was  relatively  large  at  this  point  and  the  associated  time  to 
decay  to  approximate  steady-state  conditions  would  invariably  exceed  the 
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allocated  computational  resources  on  the  CDC  7600,  the  decision  was  made 
to  both  verify  and  demonstrate  that  convergence  could  be  accelerated  by 
actively  modifying  the  governing  parameters  affecting  both  stability  and 
convergence.  The  donor-cell  and  over-rolaxation  parameters  as  well  as 
the  temporal  increment  were  reduced  to  0.9,  1.35,  and  0.0001  seconds 
respectively.  The  corresponding  step  change  in  the  incremental  value  of 
the  domain  kinetic  energy  depicted  in  Figure  17  resulted. 

IV. 2.2  Numerical  Solution  :  Results 

The  resulta  of  both  numerical  simulations  of  the  high  temperature, 
low  molecular  weight  exhaust  gases  injected  into  and  interacting  with 
the  high-subsonic  ambient  crossflow  satisfied  all  espoused  research 
objectives,  namely,  the  methodology  not  only  provided  the  means  to 
determine  the  trajectory  of  the  exhaust  plume,  but,  also  to  determine 
the  location  and  extent  of  any  recirculation  region  as  well  as  any 
region  with  possibly  large  rates  of  heat  transfer.  The  coarse 
simulation  verified  the  mathematical  model  and  the  numerical  algorithm; 
the  refined  mesh  simulation  provided  the  enhanced  data  which  form  the 
basis  of  the  results  presented  herein.  A  welcome,  but  totally 
unanticipated  result  of  both  simulations,  was  the  mathematical 
verification  of  a  flow  interaction  phenomenon  which  was  experimentally 
observed  and  also  equally  unanticipated.  Thus  the  simplified 
methodology  provided  a  mathematical  tool  of  sufficient  accuracy  which 
was  capable  of  duplicating  phyaical  phenomena  and  thereby  provided  the 
means  to  study  the  nature  of  violent  injection  problems. 


Trajectory  Analysis 


As  the  low  density  jet  impacted  the  ambient  crossflow,  only 
moderate  penetration  was  achieved,  namely  the  jet  centerline  remained 
within  one  D  of  the  injection  plane  for  the  range  of  interest  of 
distances  aft  of  the  jet  orifice*  Consequently,  the  mixing  processes 
associated  with  the  jet-crossflow  interaction  are  confined  to  a  region 
within  close  proximity  to  the.  injection  surface.  To  understand  the 
physical  nature  of  these  processes,  contour  plots  of  the  gas  mixture 
density,  temperature,  exhaust  gas  species  mass  fraction,  and  molecular 
weight  were  utilized.  From  these  isopleth  presentations,  trajectory 
data  were  easily  obtained  from  which  comparisons  to  existing  empirical 
formulae  were  made. 

Figure  18  is  a  density  contour  map  of  the  symmetry  plane  from  which 
it  is  observed  that  as  the  jet  enters  the  crossflow,  the  crossflow 
impacts  the  windward  side  of  the  jet  and  is  compressed  ahead  of  the  jet 
orifice  by  the  "quasi-solid"  appearance  of  the  jet  to  the  crossflow. 
Because  of  this  flow  stagnation,  low  density  gas  is  permitted  to  diffuse 
upstream  of  the  jet  within  the  boundary  layer  which  is  evidenced  by  the 
contours  forming  a  "lip-like"  region  ahead  of  the  jet  orifice.  Also 
from  the  figure,  it  is  apparent  that  the  low  density  gas  mixture  is 
contained  within  a  stream  tube,  the  upper  surface  of  which  is  less  than 
0.5  D  from  the  injection  plane. 

Similar  results  are  observed  in  Figure  19  in  which  contours  of 
constant  temperature  in  degrees  Kelvin  are  depicted.  High  temperature 
(low  density)  gases  are  constrained  to  a  central  stream  tube  and  the  jet 
plume  cools  rapidly  to  essentially  freestream  values  in  approximately 
2.0  D  on  the  windward  side.  Between  the  stream  tube  and  the  injection 
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surface,  a  region  of  high  temperature  gas  exists  which  is  characterized 
by  large  temperature  gradients.  Ahead  of  the  jet,  isotherms 
corresponding  to  the  density  contours  of  the  previous  figure  show  large 
streamwise  gradients  which  are  characteristic  of  a  compression  zone. 
Thermal  diffusion  in  the  boundary  layer  occurs  ahead  of  the  jet  because 
of  the  essentially  stagnated  flow  within  this  region.  Both  Figures  18 
and  19  show  that  the  region  immediately  aft  of  the  jet  is  a  region  of 
high  temperature,  low  density  gas.  It  is  revealing  to  note  that  regions 
in  the  isotherm  presentation  on  the  windward  side  of  the  jet  centerline 
exist  at  temperatures  greater  than  that  associated  with  the  gaseous 
injectant.  The  local  high  temperature  zones  are  caused  by  conversion  of 
kinetic  to  thermal  energy  within  the  mixing  region.  Finally,  there 
appear  to  be  at  least  three  regions  of  possibly  large  rates  of  heat 
transfer  to  the  injection  plane,  namely,  the  large  gradient  regions 
ahead  and  aft  of  the  jet  orifice  as  well  as  a  large  region  aft  of  the 
orifice  which  extends  approximately  4  D  in  the  streamwise  direction. 

Because  the  exhaust  species  are  corrosive  by  virtue  of  the  fact 
that  the  mixture  contains  halogen  compounds,  it  was  imperative  to 
ascertain  the  trajectory  of  the  jet  as  well  as  the  composition  of  the 
near  field  gases  within  the  high  temperature  gradient  regions.  Figure 
20  shows  that  the  high  temperature  gradient  region  within  the  jet  wake 
is  composed  predominately  of  air.  In  a  distance  of  only  0.25  1)  aft  of 
the  jet  orifice,  the  0.05  mass  fraction  contour  for  the  exhaust  gas 
shows  increased  separation  from  the  surface  as  the  streamwise  distance 
is  increased-.  Ahead  of  the  jet,  the  species  mass  fraction  contours 
replicate  the  upstream  diffusion  evidenced  in  both  Figures  18  and  19. 


However,  unlike  the  data  displayed  in  these  figures.  Figure  20  9hows 
that  the  gaseous  injectant  is  confined  to  a  region  even  closer  to  the 
injection  surface  than  indicated  by  the  density  and  temperature 
isopleths.  This  behavior  is  indicative  of  the  fact  thac  molecular 
diffusion  although  affecting  the  jet-crossflow  interaction  does  not  play 
a  role  which  is  as  significant  as  either  the  thermal  diffusion  or 
convective  processes.  Also  aoparent  from  this  f^cure  is  the 
recirculation  region  comprised  of  a  mixture  of  exhaust  gas  and  entrained 
air  on  the  leeward  side  of  the  jet  (wake  region),  the  maximum  exhaust 
gas  mass  fraction  of  which  is  approximately  0.35. 

Figure  21  depicts  contour  data  for  the  molecular  weight  of  the  gas 
mixture.  It  is  obvious  that  the  low  density  (low  molecular  weight) 
exhaust  gas  is  indeed  constrained  to  a  small  diameter  stream  tube  which 
undergoes  rapid  cooling.  The  high  temperature  region  aft  of  the  jet  and 
along  the  injection  surface  is  comprised  of  heavy  molecular  weight  gas 
(air)  in  which  large  gradients  of  molecular  weight,  in  addition  to  those 
of  temperature,  exist.  Thus,  it  is  apparent  that  molecular  diffusion 
occurs  to  a  greater  extent  within  the  wake  recirculatory  region  than  it 
does  on  the  windward  side  of  the  exhaust  plume  where  convection  is 
dominant.  Consequently,  the  nature  of  the  low  Q  (0.15)  injection  is 
predominately  controlled  by  convective  processes  with  minimal 
contributory  affects  from  molecular  diffusion  upon  the  trajectory  of  the 
exhaust  gas  jet  plume. 

The  jet  trajectories  determined  from  the  density,  temperature,  gae 
species  mass  fraction,  and  molecular  weight  coincide  as  indicated  in 
Figure  22.  In  fac*",  the  composite  trajectory  curve  appeared  to  vary  in 


74 


accordance  with  experimental  observations  and  experimentally  derived 
empirical  expresions,  that  is,  Z/D  is  proportional  to  X/D  raised  to  the 
one-third  power;  similarly,  the  slope  varies  with  X/D  raised  to  the 
negative  two-thirds  power.  Figure  23  depicts  isopleths  of  constant  Mach 
number  for  the  symmetry  plane.  These  contours  show  the  same  slope  as 
the  superimposed  trajectory  curve  determined  from  Figures  18  through  21. 
It  is  interesting  to  note  that  the  jet  centerline  curve  closely 
approximated  the  Mach  0.24  contour,  a  value  which  is  identical  to  that, 
of  the  jet  injection  Mach  number. 

Flowfield  Analysis 

In  addition  to  the  determination  of  the  jet  plume  trajectory  and 
the  regions  of  possibly  high  heat  transfer  rates  from  the  analyses  of 
the  contour  data  presented  in  Figures  18  through  23,  further  information 
was  obtained  from  analyses  of  the  flowfield  itself  and  its  Inherent 
properties.  Figure  24a  depicts  the  velocity  field  of  the  symmetry  plane 
of  the  jet  from  which  the  extent  of  the  wake  recirculation  zone  is 
easily  obtained.  This  region  which  extends  from  the  aft  portion  of  the 
jet  orifice  to  a  position  1.25  D  further  aft  along  the  injection  surface 
does  not  extend  normally  beyond  D/4  from  the  injection  surface. 
Consequently,  coupled  with  the  large  temperature  gradient  prevalent  in 
this  region,  the  recirculation  zone  further  augments  the  heat  transfer 
from  the  gas  to  the  surface.  At  the  position  3.2  D  from  the  center  of 
the  jet,  a  wall  vortex  (recirculation  zone)  formed  and  extended  to  a 
point  on  the  surface  4.5  D  from  the  jet  center.  As  observed  from 
Figures  24b  through  24c  which  are  adjacent  y-planar  cuts  of  the 
computational  domain  situated  one  and  two  cells  respectively  from  the 
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symmetry  plane  (Figure  24a),  the  lateral  span  of  the  secondary 
recirculation  zone  does  not  exceed  0.3  D.  From  the  knowledge  of  fluid 
dynamic  flows  of  single  species  perfect  gases,  a  recirculation  zone  only 
occurs  in  the  presence  of  an  adverse  pressure  gradient.  To  verify  that 
the  mathematically  predicted  wall  vortex  was  not  the  result  of  spurious 
computations  due  to  error  uncertainties  of  the  algorithm,  the  pressure, 
temperature,  molecular  weight,  and  ratio  of  specific  heats  in  the 
symmetry  plane  along  the  surface  of  the  injection  plane  were  plotted  in 
Figure  25.  As  evidenced  by  this  figure,  the  pressure  increased  upstream 
of  the  jet  and  subsequently  fell  rapidly  aft  of  the  jet  orifice  from 
which  point  the  flowfield  attempted  to  recover  from  the  pressure  defect. 
Approximately  2.2  D  aft  of  the  jet  center,  an  adverse  pressure  gradient 
formed  thereby  providing  the  necessary  requirement  for  the  development 
of  any  recirculation  zone.  The  temperature  along  the  adiabatic  wall 
increased  impulsively  from  the  freestream  value  ahead  of  the  jet  orifice 
to  a  local  maximum  at  the  jet  exit.  The  temperature  then  dropped 
rapidly  to  50  degrees  Kelvin  above  ambient  at  a  point  2.2  D  from  the  jet 
center  which  coincided  with  the  strearawise  location  of  the  adverse 
pressure  gradient.  However,  temperatures  greater  than  the  melting 
temperature  of  aluminum  alloy  (923  °K)  employed  in  aircraft  construction 
persisted  0.2  D  aft  of  the  diffuser  exit.  The  molecular  weight 
variation  exhibited  the  same  behavior  with  the  jet-crossflow  mixing  aft 
of  the  jet  along  the  injection  surface  being  essentially  completed 
within  1.4  D  of  the  jet  center  after  which  the  gas  mixture  tended  to 
recover  to  ambient  values.  The  specific  heat  ratio  variation  mirrored 
the  molecular  weight  variation  and  recovered  to  a  value  representative 
of  pure  air  by  1.4  D  aft  of  the  jet  center. 
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The  presence  of  this  wall  vortex  which  had  been  observed  to  a 
smaller  scale  in  the  less  accurate  coarse  grid  simulation  was  totally 
unanticipated.  A  similar  occurrence  had  been  observed  by  Zakanycz  (Ref 
50)  in  experimental  binary  gas  mixing  research  for  a  coflowing  helium 
jet  injected  into  air.  The  presence  of  the  experimentally  occurring 
wall  vortex  was  observed  only  after  numerical  correlation  data  for  the 
measured  turbulence  intensities  exhibited  severe  discrepancies.  As  will 
be  shown  in  subsequent  hot  air  and  helium  injectant  simulations,  this 
phenomenon  appears  to  occur  only  when  gases  of  two  vastly  different 
molecular  weights  are  mixed.  Recurrence  of  this  wall  vortex  formation 
was  observed  for  the  hot  helium  jet  but  not  for  the  hot  air  jet. 
Whether  a  similar  phenomenon  occurs  for  a  more  heavy  molecular  weight 
gas  injected  into  a  lighter  molecular  weight  gas  is  pure  conjecture  at 
this  time  as  is  the  range  of  molecular  weights  for  which  this  phenomenon 
is  a  function.  However,  the  presence  of  thi3  wall  vortex  embedded 
between  the  wall  and  the  high  temperature,  large  gradient  region  in  the 
jet  wake  could  add  an  appreciable  increment  to  the  heat  transfer  in  this 
region  computed  in  any  simplified  analysis  which  did  not  account  for  its 
presence.  If  the  aft  portion  of  the  injection  surface  was  modified  to 
simulate  the  curvature  of  the  fuselage  of  the  wide-bodied  air  vehicle, 
it  is  doubtful  that  the  wall  vortex  would  disappear  because  increased 
curvature  introduces  additional  buoyancy  effects  into  the  flowfield 
thereby  almost  positively  assuring  that  some  vortex  region  would  exist 
despite  the  pressure  relief  afforded  by  the  increased  curvature. 

Figure  26  presents  pressure  coefficient  data  for  the  symmetry 
plane.  From  observation  of  the  data,  there  exist  two  main  pressure 
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centers.  In  the  compression  zone  ahead  of  the  jet,  a  large  positive 
pressure  cell  exists  approximately  D/2  ahead  of  the  jet  and  D/2  down 
from  the  injection  surface.  A  smaller  negative  pressure  cell  exists  in 
the  wake  region  at  the  aft  edge  of  the  jet  orifice  (D/2)  and  apparently 
centered  on  or  very  close  to  the  injection  surface  itself.  This  low 
pressure  region  entrains  high  temperature  waste  gases  thereby  creating 
the  high  heat  transfer  which  was  identified  earlier  from  trajectory 
contour  data. 

iv. 3  Simplified  Analyses 

Because  of  the  exhorbitant  amount  of  computational  resources 
required  to  solve  the  three-dimensional  Navier-Stokes  for  real  gas 
mixtures  with  varying  transport  and  thermodynamic  properties,  two 
additional  simulations  were  made  to  determine  whether  simplified 
analyses  could  be  utilized  to  obtain  usefull  engineering  data  which 
would,  in  fact,  replicate  the  data  obtained  from  more  costly  gas  mixture 
analyses.  These  simplified  analyses  included  a  hot  helium  jet  and  a  hot 
air  jet  constrained  by  an  injection  temperature  of  1500  °K,  an  injection 
Mach  number  of  0.24,  and  an  injection  to  freestream  dynamic  pressure 
ratio  of  0.15.  Additionally,  it  was  further  anticipated  that  the 
perfect  gas  simulation  for  hot  air  injected  into  the  cold  ambient 
crossflow  would  provide  some  means  of  correlating  real  gas  mixture 
injection  with  chat  of  the  more  simply  modeled  hot  air  injection. 

IV. 3.1  Helium  Jet 

The  simulation  of  an  aspect  ratio  1.75  rectangular  jet  of  pure 
helium  at  a  temperature  of  1500  °K  being  injected  into  a  Mach  0.7 


78 


crossflow  at  an  aircraft  pressure  altitude  of  11  kilometers  was 
accomplished  by  dividing  the  computational  domain  into  29  x  13  x  13 
cells  in  the  x,  y,  and  z  spatial  coordinate  directions  respectively  such 
that  the  CDC  Cyber  175  at  Wright-Patterson  AFB  was  capable  of  handling 
the  resultant  core  requirements.  A  K  value  of  1.20  was  employed  to 
expand  the  geometrically  varyirg  z-dimension  of  each  cell  from  the 
injection  surface.  The  jet  orifice  was  modeled  by  4  x  3  cells  in  the  x 
and  y  directions  in  the  injection  plane  respectively.  The  numerical 
solution  was  terminated  after  an  acceptable  level  of  convergence  had 
been  obtained  as  is  evidenced  by  Figure  27.  The  typical  approach  to 
steady-state  was  observed  and  computations  were  terminated  after  t=0.34 
seconds  which  corresponds  to  approximately  2.5  characteristic  times. 

The  flcwfield  data  of  the  symmetry  plane  are  presented  in  Figure 
28.  The  recirculation  zone  in  the  wake  region  appears  to  be  more 
extensive  than  that  predicted  by  the  simulation  utilizing  the  more 
refined  mesh.  A  minute  indication  of  the  wall,  or  secondary,  vortex 
along  the  injection  plane  at  approximately  4.0  D  from  the  jet  center  is 
also  indicated.  However,  the  coarse  grid  simulation  for  helium  fails  to 
define  the  extent  of  this  recirculation  zone. 

The  density  isopleths  in  Figure  29  depict  the  compression  zone 
ahead  of  the  jet  orifice  as  well  as  the  region  characterized  by  high 
temperatures  and  large  temperature  and  molecular  weight  gradients  in  the 
wah t  near  the  injection  surface.  The  jet  trajectory  obtained  from  this 
figure  is  plotted  in  Figure  22.  It  is  evident  that  the  hot  helium 
Injectant  trajectory  curve  replicates  the  Z/D  variation  with  X/D  to  the 
one-third  power;  however,  the  error  induced  by  the  coarse  mesh  helium 
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simulation  resulted  in  the  trajectory  being  situated  above  the 
corresponding  curve  derived  from  the  refined  mesh  ga9  mixture 
simulation.  As  the  X/D  value  increased,  the  agreement  between  these  two 
curves  increased. 

IV. 3. 2  Hot  Air  Jet 

The  simulation  of  an  aspect  ratio  1.75  rectangular  air  jet  at  1500 
°K  being  injected  into  an  ambient  crossflow  at  0.7  Mach  number  at  an 
aircraft  pressure  altitude  of  11  kilometers  was  similar.ly  accomplished 
by  dividing  the.  computational  domain  into  29  x  13  x  17  cells  in  the  x, 
y,  and  z  spatial  coordinate  directions  respectively  such  that  the 
simulation  could  be  conducted  on  the  CDC  Cyber  175  at  Wright-Patterson 
AFB.  A  K  factor  of  1.20  was  also  used  to  expand  the  z-dimension  cell 
length  for  every  cell  measured  from  the  injection  surface. 
Additionally,  the  jet  orifice  was  modeled  by  4  x  3  cells  in  the  x  and  y 
directions  respectively  embedded  in  the  injection  surface.  The 
numerical  solution  was  terminated  after  t=0.24  seconds  despite  the  fact 
that  the  solution  had  exhibited  an  acceptable  level  of  convergence  long 
before  (Figure  30)  and  no  numerical  instability  appeared  probable.  The 
termination  time  in  this  case  was  still  less  than  two  characteristic 
times. 

As  is  evidenced  by  the  flowfield  data  of  the  symmetry  plane  in 
Figure  31,  only  the  recirculation  zone  within  the  wake  was  discerned 
from  the  numerical  simulation.  As  alluded  to  in  Section  IV. 2. 2,  the 
appearance  of  the  wall  vortex  did  not  recur  which  reinforced  by  the 
behavior  of  the  flowfield  tends  to  confirm  the  hypothesis  of  the  vortex 
being  driven  by  a  species  diffusive  flux  caused  by  gradients  in 
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molecular  weight  or  species  concentrations  in  the  presence  of  a  "no- 
slip”  surface. 

The  density  contour  plot  in  Figure  32  depicts  a  similar  compression 
zone  ahead  of  the  jet  orifice  and  the  high  temperature,  large  gradient 
layer  in  the  wake  and  adjacent  to  the  injection  surface.  The  trajectory 
curve  obtained  from  this  figure  was  also  plotted  in  Figure  22.  However 
in  this  case,  the  trajectory  lies  below  that  of  the  refined  grid  gas 
mixture  simulation.  It  too  demonstrates  the  same  general  behavior  as 
the  other  trajectories,  nanely,  the  penetration  depth,  Z/D,  variation 
with  X/D  to  the  one- third  power.  Golubev  and  Klimkin  (Ref  51)  have 
shown  that  for  jets  with  equal  initial  momenta,  "the  rate  at  which  the 
apparent  additional  mass  increases  is  higher  in  jets  with  low  density 
than  in  jets  with  high  density."  Consequently,  the  lower  trajectory  for 
the  hot  air  is  thus  easily  explained  by  the  constraint  of  constant  Q 
(0.15)  which  is  the  ratio  of  jet  to  freestream  momentum  fluxes  while 
also  subject  to  the  additional  constraint  of  equal  jet  injection 
temperatures.  The  density  of  air  which  is  greater  than  that  of  helium 
requires  that  velocity  squared  terra  of  air  be  less  than  the 
corresponding  term  for  helium  to  maintain  the  constancy  of  the  density- 
velocity  squared  product  (Q) .  The  resulting  exit  velocity  of  the  air 
jet  (176  m/sec)  was  therefore  approximately  one-third  that  of  the  helium 
jet  (547  m/sec)  at  the  diffuser  exit  for  a  jet  exit  Mach  number  equal  to 
0.24.  Because  of  the  decreased  injection  velocity,  less  penetration  of 
the  crossflow  was  achieved.  Further  attenuation  of  the  jet  penetration 
was  caused  by  the  entrainment  of  the  ambient  crossflow. 
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During  the  course  of  the  research,  it  was  observed  that  the  je_- 
crossflow  interactions  of  both  low  and  high  molecular  weight  gas 
injectanta  with  the  ambient  crossflow  exhibited  essentially  the  same 
convection  dominated  behavior  and  contained  essentially  the  same 
physical  phenomena.  The  main  difference,  excluding  obvious 
concentration  induced  effects,  was  the  degree  to  which  these 
characteristics  were  observed.  Because  of  the  large  differences  in 
molecular  weights  of  the  injected  gases,  it  was  necessary  to  employ  some 
correction  factor  to  account  for  this  variation  if  the  trajectories  were 
to  collapse  into  a  single  curve. 

In  typical  heat  and  mass  transfer  problems  involving  gas  injectants 
not  of  the  same  species  as  that  of  the  dominant  gas  medium,  the  ratio  of 
the  mass  transfer  conductance  to  the  corresponding  value  with  "no 
blowing"  is  corrected  by  the  following  factor  which  is  dependent  upon 
the  flow  characteristics  (Ref  10,  pages  322  and  326): 


G  = 


freestream  molecular  weight 
considered  phase  surface  molecular  weight 


(100) 


where  s  is  either  2/3  or  0.4  for  laminar  and  turbulent  flow, 
respectively.  Though  no  logical  relation  exists  to  extend  Eq  (100)  for 
use  at  high  blowing  rates,  its  use  as  an  empirical  correction  factor 
fortuitiously  does  account  for  molecular  weight  variations.  In  the 
present  research  of  turbulent  rectangular  jet  injection,  the  value  of 
the  exponent  s  in  Eq  (100)  was  assumed  to  be  0.4  and  the  molecular 
weight  of  the  considered  phase  was  assumed  to  be  identical  to  that  of 
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the  injectant.  The  results  of  applying  this  factor  to  both  simulation 
and  empirically  derived  jet  trajectories  are  depicted  in  Figure  22.  The 
correction  factor  modification  clearly  demonstrated  the  proper  behavior; 
with  additional  analyses  and  experimentation,  enhanced  results  may  be 
obtained  by  more  accurately  describing  the  exponent.  Of  all  the 
trajectory  expressions  analyzed  during  this  research,  only  that  of 
Callaghan  and  Ruggeri,  Eq  (8),  shows  some  promise  of  being  used  in 
conjunction  with  the  molecular  weight  correction,  G,  given  by  Eq  (100) 
for  interacting  gases  with  different  molecular  weights  at  elevated 
temperatures  to  predict  the  jet  plume  centerline.  In  fact,  when  applied 
to  the  conditions  imposed  on  the  injected  gas  mixture  and  crossflow  and 
utilizing  either  the  streamwise  dimension,  D,  or  the  equivalent 
diameter,  De,  of  the  jet  orifice  in  the  expression,  the  resultant  curve 
varied  less  than  10  and  5  per  cent  respectively  from  that  derived  from 
the  refined  mesh  numerical  simulation. 

IV . 5  Turbulence  Diffusivity:  Heat  Transfer  Effects 

Additional  numerical  simulations  were  conducted  to  ascertain  the 
effects  of  turbulence  diffusivity  upon  the  heat  transfer  to  the 
injection  surface  (aircraft  fuselage).  These  effects  were  essentially 
bracketed  by  simulations  made  with  both  molecular  and  turbulent  values 
of  the  diffusivities  for  air  and  helium  gaseous  injectants  at  1500  CK 
further  constrained  to  the  same  Q  ratio,  Mach  numbers  associated  with 
the  jet  and  freestream,  and  operating  altitude  of  the  previous  studies. 
These  sensitivity  analyses  were  accomplished  on  the  CDC  Cyber  750 
(replacement  for  the  CDC  Cyber  175)  at  Wright-Pattersoa  AFB,  Ohio, 
utilizing  the  identical  coarse  meshes  described  earlier.  The  "no-slip" 
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injection  surface  was  defined  to  be  adiabatic,  that  is,  no  heat  transfer 
to  the  surface  from  the  gas  was  permitted;  however,  the  adiabatic  wail 
temperature  does  provide  an  accurate  indication  of  the  heat  transfer  to 
the  surface  as  it  has  been  sufficiently  demonstrated  in  the  literature 
pertaining  to  film  cooling.  In  each  of  the  ensuing  simulations,  the 
turbulence  model  utilized  was  that  given  by  Eq  (34).  The  minimal 
diffusivity  runs  employed  the  effective  viscosity  in  concert  with  just 
the  molecular  (laminar)  values  of  thermal  conductivity  and  mass  transfer 
binary  diffusion  coefficient  (helium  runs,  only).  The  turbulent 
diffusivity  runs  incorporated  the  effective  values  of  viscosity,  thermal 
conductivity,  and  binary  diffusion  coefficient. 

The  adiabatic  wall  temperatures  for  air  with  minimal  and  turbulent 
diffusivities  are  presented  in  contour  map  form  in  Figures  33  and  34 
respectively  from  which  the  nature  of  the  heat  transfer  to  the  surface 
can  be  directly  inferred.  When  just  the  molecular  value  of  thermal 
conductivity  was  used  (Figure  33)  only  a  near  field  thermal  diffusive 
effect  exists  such  that  the  wall  temperature  recovers  to  essentially 
frees tream  values  within  D/4  of  the  diffuser  exit.  However,  when  the 
full  turbulent  dif fusivities  were  employed  (Figure  34),  the  thermal 
diffusion  was  significantly  enhanced.  High  temperature  (>550  °K) 
regions  exist  approximately  D/4  upstream  of  the  diffuser  exit,  D/2 
laterally  from  the  diffuser  exit,  and  3D  aft  of  the  diffuser  exit  along 
the  center  line.  This  drastically  large  region  of  elevated  wall 
temperature  poses  severe  thermal  protection  problems  which  must  be 
solved. 

The  corresponding  adiabatic  wall  temperature  presentations  for 


helium  with  both  minimal  and  fully  turbulent  dif fusivities  are  presented 
in  Figures  35  and  36,  respectively.  The  solution  for  the  minimal 
diffusivity  run  was  still  developing  when  the  simulation  was  terminated. 
However,  the  extent  of  the  high  temperature  region  is  far  greater  than 
the  associated  region  for  air  (Figure  33) .  When  the  full  set  of 
dif fusivities  was  employed  as  evidenced  in  Figure  36,  the  increased 
thermal  diffusion  is  even  more  drastic  than  that  produced  by  the 
turbulent  diffusivity  run  for  air.  Not  only  does  the  high  temperature 
region  encompass  more  of  the  injection  surface,  but  also,  it  is  evident 
that  an  exceptionally  large  region  exists  in  which  the  local  temperature 
(adiabatic  wall)  is  in  excess  of  the  melting  temperature  of  aluminum 
alloys  employed  in  aircraft  construction.  This  melting  region  is 
contained  approximately  within  the  subregion  D/4  upstream  of  the 
diffuser  exit,  D/2  laterally  from  the  diffuser  exit,  and  2.5  D  aft  of 
the  diffuser  exit.  As  is  evidenced  by  Figure  36,  surrounding  this 
melting  region  a  far  larger  region  exists  in  which  the  aluminum  alloys 
would  lose  significant  percentages  of  the  available  structural  strength 
and  consequently  would  lead  to  structural  deformation.  Because  the  AFWL 
gas  mixtures  are  comprised  of  approximately  86  per  cent  helium,  the 
temperature  distribution  on  the  injection  surface  (fuselage)  and  the 
associated  heat  transfer  to  it  should  virtually  replicate  the 
corresponding  behavior  and  characteristics  associated  with  the  pure 
helium  jet  injection. 

Based  upon  these  diffusivity  sensitivity  studies,  the  heat  transfer 
mechanism  to  the  injection  surface  can  be  divided  into  two  coupled 
facets.  First,  because  of  the  enhanced  dif fusivities  provided  by  the 
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turbulent  transport,  significant  thermal  diffusion  occurs  between  the 
jet  plume  and  the  injection  surface.  However,  with  an  increase  in 
turbulent  diffusivity,  the  second  facet  of  the  heat  transfer  mechanism 
comes  into  play.  The  high  turbulent  diffusivities  associated  with  the 
jet  plume  cause  the  plume  to  noth  expand  and  cool  thereby  directly 
reducing  the  source  of  energy  for  heat  transfer  to  the  injection 
surface.  Consequently,  the  wall  temperature  then  decreases  to  somewhat 
lower  levels  because  of  both  the  enhanced  turbulent  diffusive  efficiency 
and  the  reduction  of  available  energy  within  the  jet  plume. 


IV. 6  Summary 


The  hot  (1500  °K)  exhaust  gases  from  an  airborne  chemical  laser 
being  ejected  downwards  from  an  aspect  ratio  1.75  rectangular  jet 
orifice  aligned  with  the  major  axis  parallel  to  the  streamwise  Mach  0.7 
ambient  crossflow  was  successfully  simulated.  The  data  acquired  from 
the  numerical  simulation  described  the  jet  plume  trajectory,  the  extent 
of  the  recirculation  zone  in  the  wake  of  the  jet-crossflow  interaction, 
and  the  regions  of  high  possible  rates  of  heat  transfer.  Additionally, 
the  numerical  method  demonstrated  the  existence  of  a  secondary  wall 
vortex  (recirculation  zone)  which  formed  after  the  jet  plume  was 
sufficiently  deflected  to  be  essentially  aligned  with  the  ambient 
crossflow  (coflowing).  The  formation  of  a  wall  vortex  was 
experimentally  observed  by  Zakanycz  (Ref  50)  in  his  experimental 
investigation  of  binary  ga3  mixing.  Nowhere  else  In  the  available 
literature  has  this  phenomenon  been  mentioned  which  apparently  only 
occurs  for  low  molecular  weight  gases  interacting  with  a  higher 
molecular  weight  coflowing  gaseous  medium  (air,  for  example)  in  the 
presence  of  a  constraining  ''nc-slip"  surface. 
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Simplified  analyses  were  also  conducted  to  determine  if  the 
essential  flow  phenomena  were  captured  mathematically  utilizing  a  simple 
binary  gas  interaction,  that  is,  hot  helium  injected  into  croesf lowing 
ambient  air,  or  even  more  simply,  just  utilizing  a  hot  air  jet  being 
injected  into  the  ambient  crossflow.  The  essential  flow  features  of  the 
exhaust  gas  mixture  simulation  were  replicated  by  the  helium  jet 
simulation  to  some  extent;  however,  the  accuracy  of  the  coarse  grid 
simulation  was  not  to  the  same  level  of  accuracy  provided  by  the  refined 
mesh  gas  mixture  simulation.  Additionally,  the  wall  vortex  formation 
which  did  become  apparent  in  the  felium  injectant  simulation  failed  to 
recur  for  the  hot  air  jet-crossflow  interaction  problem  thereby  adding 
credence  to  the  hypothesis  that  diffusive  fluxes  induced  the  formation 
of  this  secondary  wall  vortex. 

Not  only  were  useful  engineering  data  obtained  from  the  simplified 
numerical  algorithm,  but,  also  the  utility  of  the  methodology  as  a 
learning  and  prognostic  tool  to  study  the  effects  of  violent  jet- 
crossflow  interactions  was  demonstrated.  The  flowfields  for  both  real 
gases  and  hot  air  showed  essentially  the  same  basic  structure  and 
characteristics  because  of  the  dominance  of  the  convective  processes. 
However,  the  effects  of  species  concentration  gradients  which  apparently 
induced  a  secondary  recirculation  zone  to  form  could  not  be  duplicated 
by  utilizing  only  one  gas  in  the  jet-crossflow  interaction  problem. 
Sufficient  evidence  was  obtained  to  warrant  the  simplification  of  a  gas 
mixture  injection  problem  by  replacing  the  mixture  with  a  representative 
gas  such  that  an  effective  binary  gas  mixture  in  the  flowfield  would 
result  with  the  attendant  savings  in  computational  resources;  if 
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chemical  reactions  were  permitted  to  occur  within  the  simulation,  this 
simplification  could  not  obviously  be  made.  The  coarse  grid  utilized  in 
the  simplified  analyses  resulted  in  numerical  errors  which  were 
unacceptable  for  quantitative  descriptions  of  the  observed  flow 
phenomena  but  which  were  acceptable  for  qualitative  investigations.  For 
results  of  engineering  utility,  a  finer  subdivision  of  the  computational 
domain  is  required  and  a  geometric  stretch  factor,  K,  not  exceeding  1.09 
should  be  used  because  of  the  inherent  inaccuracy  of  the  numerical 
algorithm  when  a  non-uniform  mesh  is  employed. 

The  application  of  a  molecular  weight  correcton  factor  given  by  Eq 
(100)  to  the  simplified  analyses  demonstrated  that  trajectory  curves 
could  be  synthesized  from  one  another  for  identical  jet  and  crossflow 
contraints.  Use  of  this  factor  in  conjunction  with  the  empirical 
trajectory  expression  of  Callaghan  and  Ruggeri  resulted  in  less  than  5 
per  cent  deviation  from  the  trajectory  derived  from  the  refined  mesh 
numerical  simulation  of  the  gas  mixture  (Table  I)  injectant. 

Finally,  the  effects  of  the  turbulent  dif fusivities  upon  the  heat 
transfer  to  the  injection  surface  were  determined  for  both  air  and 
helium  gas  injection  problems.  Increased  levels  of  heat  transfer  to  the 
surface  resulted  when  increased  levels  of  dif fu3ivities  were  utilized 
(as  indicated  bv-  the  adiabatic  wall  temperatures).  Regions  of  severe 
heating  were  found  to  exist  for  both  high  temperature  air  and  helium 
injections  at  low  values  of  Q.  The  heating  was  especially  pronounced 
for  the  helium  case  where  large  surface  regions  were  discerned  to  exist 
in  which  the  injection  surface  adiabatic  wall  temperatures  exceeded  both 
the  temperature  associated  with  structural  deformation  and  the  melting 
temperature  of  aluminum  alloys  employed  in  aircraft  construction. 
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V.  Conclusions  and  Recommendations 


The  high  temperature  exhaust  gases  from  an  airborne  chemical  laser 
ejected  normally  downward  from  an  aspect  ratio  1.75  rectangular  diffuser 
orifice  mounted  in  the  lower  surface  of  the  fuselage  of  a  wide-bodied 
aircraft  into  the  Mach  0.7  ambient  crossflow  at  an  aircraft  pressure 
altitude  of  11  kilometers  was  successfully  numerically  simulated.  A 
simplified  finite  difference  donor-cell  ICE  algorithm  (Ref  44)  was 
modified  to  solve  the  complete  set  of  time  dependent,  compressible, 
three-dimensional  Navier-Stokes  equations  in  addition  to  a  species 
conservation  equation  while  employing  variable  thermodynamic  and 
transport  properties.  Turbulence  closure  was  achieved  through  the  use 
of  a  locally  varying  velocity  defect  eddy  viscosity  model.  In  the 
mathematical  model  of  the  physical  problem,  the  gas  injectant  was 
proscribed  from  reacting  with  the  ambient  crossflow.  (C-l)  The  numerical 
methodology  exhibited  the  capability  to  solve  violent  gaseous  injection 
problems  such  as  the  normally  ejected  jet  interacting  with  an  ambient 
crossflow  at  both  high,  and  more  importantly,  low  ratios  of  jet  to 
crossflow  momentum  flux  (16.0  and  0.15). 

Constrained  by  the  same  injection  Mach  number,  Q  value,  and 
temperature,  (C-2)  the  basic  structure  and  characteristics  of  the  jet- 
crossflow  interaction  with  the  Mach  0.7  crosssflow  at  a  simulated 
pressure  altitude  of  11  kilometers  (tropopause)  were  essentially  the 
same  regardless  of  the  jet  injection  species.  (C-3)  Convective  processes 
were  seen  to  dominate  the  jet-crossflow  interaction  even  when  helium  was 
the  injectant.  (C-4)  Thermal  diffusion  was  seen  to  have  a  significantly 
greater  effect  on  the  resultant  gas  mixture-air  flowfield  uhan  molecular 
diffusion. 
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The  (C-5)  penetration  of  the  impacting  crossflow  by  the  fully 
turbulent  gas  injectants  when  the  Q  value  was  constrained  to  0.15  was 
inversely  proportional  to  the  densities  of  the  injectants.  Thus,  a 
helium  Jet  penetrated  the  crossflow  farther  than  the  gas  mixture 
(molecular  weight  6.6)  which  in  turn  penetrated  farther  than  the 
corresponding  air  jet.  Because  the  value  of  Q  was  constant,  the 
velocity-density  squared  product  of  the  jet  species  prior  to  injection 
was  also  fixed  thereby  requiring  the  velocity  squared  term  to  vary 
accordingly  to  satisfy  the  imposed  constraints  (Q  and  Mj). 
Consequently,  the  greater  velocity  (momentum)  defect  of  the  heavier 
molecular  weight  injectants  resulted  in  proportionately  less  crossflow 
penetration.  Additionally,  ( C— 6 )  when  the  gas  mixture  injectant  which 
was  assumed  to  be  unreactive  with  the  ambient  crossflow  was  replaced  by 
a  single  representative  gas  (for  example,  helium)  to  reduce  the 
computational  complexity,  the  algorithm  produced  results  which  were  in 
agreement  with  those  generated  by  employing  the  gas  mixture. 

As  is  evidenced  by  Figure  22,  (C-7)  utilization  of  the  molecular 


weight  correction 

factor 

described  in 

Section  IV. k 

enabled 

the 

trajectory  curve 

of  one 

species  to 

approximately 

generate 

the 

corresponding  trajectory  of  another.  (C-8)  The  empirical  relationship 
derived  by  Callaghan  and  Ruggeri,  Eq  (8),  was  employed  in  conjunction 
with  the  molecular  weight  correction  factor,  Eq  (100),  to  predict  the 
trajectory  of  the  real  gas  mixture  jet  being  ejected  into  the  ambient 
crossflow.  Excellent  agreement  was  obtained  between  the  refined  mesh 
numerical  simulation  and  the  empirically  derived  trajectories- 

Sensitivity  studies  conducted  to  determine  the  effects  of  the 
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magnitude  of  turbulent  dif fusivities  on  the  heat  transfer  to  the 
injection  surface  demonstrated  that  (C-9)  large  regions  exist  on  the 
injection  surface  in  which  temperatures  exceed  the  maximum  safe 
temperature  of  aluminum  alloys  associated  with  the  onset  of  permanent 
structural  deformation.  Severe  heating  regions  within  this  zone  have 
temperature  levels  in  excess  of  the  melting  temperature  of  aircraft 
aluminum  alloys.  Furthermore,  it  was  observed  that  (C— 10)  the  heat 
transfer  mechanism  from  the  jet  plume  to  the  injection  surface  was 
comprised  of  two  coupled  phenomena  which  lead  to  significant  levels  of 
injection  surface  heating.  As  the  turbulent  dif fusivities  were 
increased,  the  heat  transfer  to  the  surface  increased  until  a  local 
maximum  was  obtained,  after  which,  additional  increases  in  diffusivity 
resulted  in  lower  heat  transfer  rates  to  the  injection  surface  by  the 
mechanism  described  in  Section  IV. 5. 

The  following  recommendations  for  further  study  are  submitted  to 
provide  further  refinement  to  (1)  the  numerical  methodology  with  the 
attendant  increases  in  run  time  per  cycle  and  fiscal  charges,  and,  to 
(2)  the  mixing  mechanisms  involved  in  high  temperature,  low  molecular 
weight  jet  injection  into  crossflowing  gaseous  mediums  of  the  same  or  of 
a  different  species: 

(R—  1 )  Diffusion  -  add  both  the  DuFour  (’'diffusion- 
thermo")  and  the  Soret  ("thermal-diffusion")  effects  to  the 
numerical  algorithm.  Because  the  temperature  of  the  jet  and 
the  associated  interacting  flow  is  extremely  high  with 
correspondingly  large  temperature  gradients ,  the  inclusion  of 
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these  "secondary"  effects  may  have  a  pronounced  effect  on  the 
computational  results. 

(R-2)  Coordinates  -  employ  a  generalized  body  oriented 
coordinate  (geometry)  subroutine  replete  with  the  requisite 
metrics  to  permit  calculations  from  arbitrarily  defined 
bodies. 

(R-3)  Species  Generation  -  permit  the  exhausting  gas  to 
chemically  react  with  the  ambient  crossflow  within  the  wake 
region  of  the  jet-crossflow  interaction. 

(R-4)  State  -  modify  the  solution  algorithm  to  account 
for  gradients  in  concentration  (mass  fraction  and  molecular 
weight)  as  well  as  temperature.  The  existing  algorithm 
accounts  for  small  variations  in  pressure  in  the  ICE  iteration 
while  essentially  neglecting  the  variations  of  temperature  and 
concentration.  Incorporation  of  the  additional  terms  in  the 
code  would  definitely  increase  its  robustness. 

(R-5)  Radiation  -  incorporate  the  effects  of  radiative 
heat  transfer  from  the  high  temperature  boundary  layer  in  the 
wake  of  the  jet. 

(R-6)  Surface  Heat  Transfer  -  remove  the  constraint  of 
an  adiabatic  wall  to  permit  calculation  of  the  effects  of 
various  gaseous  injectants  on  the  heat  transfer  to  the  air 
vehicle  surface. 

(R-f)  Self -Optimization  -  employ  a  subroutine  in  the 
algorithm  to  automatically  compute  and  update  the  donor-cell 
parameter,  over  relaxation  factor,  and  temporal  increment 
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during  each  time  step  to  reduce  the  computational  time 
required  to  obtain  the  asymptotic  approach  to  steady-state  of 
the  numerical  solution. 

(R-8)  Molecular  Weight  -  conduct  experimental  jet- 
crossflow  binary  gas  studies  to  obtain  quantitative  data 
relating  the  trajectory  of  the  jet  plume  with  the  associated 
molecular  weights  of  the  interacting  species. 

(R-9)  Accuracy  -  employ  a  higher  order  algorithm  in 
conjunction  with  a  substantially  refined  mesh  in  the  near 
field  region  of  the  jet.  This  recommendation  would  require 
the  availability  of  a  computer  with  vastly  greater  core 
capacity  than  that  of  the  CDC  7600. 
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Domain:  Aircraft  Jet  Injection  Model 
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Figure  4.  Finite  Subdivision  of  Computational  Domain 
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Figure  6.  Computational  Domain:  Weston  and  Thames  Experimental  Model 


Figure  7.  Convergence  Criterion:  Weston  Simulation t  K=l.l,  Coarse  Mesh 


Figure  8.  Convergence  Criterion:  Weston  Simulation,  K=1.05,  Coarse  Mesh 
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Figure  9.  Convergence  Criterion:  Weston  Simulation,  K=1.05,  Refined  Mesh 
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Figure  10b.  Y-Plane  Flowfield:  Symmetry  Plane 
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Figure  12.  Symmetry  Plane  Flowfield:  Coarse  Mesh 


Figure  13a.  Pressu 
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Figure  14a.  Peake  Experiment  Pressure  Coefficient  Contour  Map,  R*-4.0 


Figure  15.  Surface  Pressure  Distribution 
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Figure  19.  Trajectory  Analysis:  Temperature  Contour  Map,  AFWL  Gas 


Figure  20.  Trajectory  Analysis:  Mass  Fraction  Contour  Map,  AFWL  Gas 


Figure  22.  Jet  Trajectory  Coraparis 
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Figure  24a.  Flowfield  Analysis:  Symmetry  Plane,  AFWL  Gas 
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Figure  24c.  Flowfield  Analysis:  Symmenry  Plane  -2,  AFWL  Gas 


Flowfield  Analysis:  Symmetry  Plane  Pressure  Coefficient  Map 


Converge 


Flowfield  Analysis:  Helium,  K— 1«20,  Coarse  Mesh 


Analysis:  Helium,  K=1.2%  Coarse  Mesh 


Figure  30.  Convergence  Criterion:  Hot  Air,  K=1.20,  Coarse  Mesh 
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Figure  31.  Flowfield  Analysis:  Hot  Air,  K=-1.20,  Coarse  Mesh 
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Figure  33.  Adiabatic  Wall  Temperature:  Hot  Air,  K=*1.20,  Coarse  Mesh,  Molecular 


Figure  34.  Adiabatic  Wall  Temperature:  Hot  Air,  K=1.20,  Coarse  Mesh,  Turbulent 
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20,  Coarse  Mesh,  Molecular 


abatic  Wall  Temperature:  Helium,  K=1.?0,  Coarse  Mesh,  Turbulent 


Table  1 


AFWL  Gaa  Mixture  (CASE  2G) 


Speciesj 

Mole  Fraction 

[Xj3 

Molecular  Weigl 

CF4 

0.01281 

160.39 

He 

0.86242 

4.0026 

°2 

0.04140 

4.028204 

DF 

Of  03852 

21.0125 

HF 

0.02562 

20.006303 

n2 

C.  01923 

2d. 016 

Average  Gas  Temperature:  1500  °K 
Average  Gas  Molecular  Weight:  6.60674 
Specific  Heat  Ratio:  1.5484 


1962  United  States  Standard  Atmosphere 


'.(he  computer  simulation  utilized  a  numerical  description  of  the 
first  two  layers  of  the  1962  standard  atmosphere.  Figure  A1  depicts  the 
temperature  variation  with  altitude  as  adopted  by  the  1962  standard 
atmosphere  (Ref  A.l). 


Figure  Al.  1962  Standard  Atmosphere 


Associated  with  the  temperature  variation  prescribpd  by  Figure  Al,  there 
were  values  of  basic  constants  utilized  by  mathematical  model  of  the 
atmosphere  which  were  also  prescribed  by  the  committee  which  derived  the 
standard  atmosphere.  The  following  list  of  constants  required  by  the 
atmospheric  model  are  as  follows: 


Universal 

Gas  Constant,  Rgflg 

1545.31 

ft-lbf/(lbm~mole-°R) 

Sea  Level 

Pressure,  Pg^ 

2116.22 

lbf/ft2 

See  Level 

Temperature,  Tg^ 

288.15 

°K 

Acceleration  of  Gravity,  g  32.174049  ft/sec^ 

Molecular  Weight  of  Air,  M  28.9644 

To  further  reduce  the  complexity  of  the  mathematical  description  of  the 
atmosphere,  non-dimensional  quantities  were  defined  to  make  the 
numerical  model  less  dependent  upon  the  myriad  of  units  which  could 
possibly  be  used  to  parameterize  the  atmosphere.  These  non-dimensional 
quantities  are  given  by  the  following: 

(D  0  -  Tambient2TSL  *  Tambla„£/288- 15 
<2>  4  '  pan,bie„t/pSL  '  pa»bl«,t/2116-22 
<3>  0  -p  ambient  SL  *  p  ambient'0-00237691 

A.l  Temperature  Ratio 

The  first  layer  of  the  atmosphere  is  defined  to  exist  between 
geopotential  altitudes  of  0  through  11  kilometers  such  that  a  uniform 
lapse  rate,  Lp  of  -6.5  °K/km  is  maintained.  The  temperature  at  any 
geopotential  altitude  between  the  endpoints  of  this  layer  is  simply 
expressed  by  the  following  equation, 

Tambient  =  TSL  "  L1  *hp  (A-1) 

where  hp  is  any  arbitrary  layer-1  altitude. 

Non-dimensionalization  of  the  above  expression  is  accomplished  by 
dividing  by  the  base  temperature,  Tg^,  that  is, 
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®  1  “  ^ambient  ”  *  L1 


1  "  K1  hp 


(A-2) 


where  *  6.8755856(10)“°  for  hp  given  in  feet. 


Noting  that  the  second  layer  is  isothermal,  that  is,  the  lapse  rate  L2 
is  identically  zero,  the  non-dimensional  temperature  ratio  between 
geopotential  altitudes  of  11  and  20  kilometers  is  a  constant  given  by. 


'2  =  216.65/288.15  =  0.75186535 


(A-3) 


The  temperature  ratio  specified  by  Eqs  (A-2)  and  (A-3)  completely  define 
the  temperature  variation  within  the  first  two  layers  of  the  standard 
atmosphere. 


A. 2  Pressure  and  Density  Ratios 

The  temperature  variation  for  Layer-1  is  prescribed  precisely  by  Eq 
(A-l).  However,  the  remaining  state  variables  required  to  describe  the 
atmosphere  still  must  be  determined >  Assuming  that  the  atmosphere 
remains  in  hydrostatic  equilibrium,  the  associated  equation  which 
relates  pressure  changes  to  corresponding  changes  in  altitude  is  given 


dp  -  -P  g  dhD 


(A-4) 


Substituting  Eqs  (A-l),  (A-2),  and  (A-4)  into  the  equation  of  state 
results  in  the  following  ordinary  first  order  differential  equation 
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g 


(A-5) 


dp 

P 


dh 


1.8  EIsl1-  Kj  hp 


where  R  -  Rgas/M* 


which  is  easily  integrated  between  sea  level  and  the  tropopause  (11  km) 
to  yield  the  following  non-dimensional  expression  for  the  pressure 
variation  within  the  first  atmospheric  layer  in  terms  of  the  non- 
dimensional  temperature  ratio, 


-  05. 2559121 


(A— 6 ) 


The  non-dimensional  pressure  variation  within  the  second  layer  of 
the  standard  atmosphere  is  obtained  circuitously  by  substituting  the 
value  of  the  temperature  of  this  isothermal  layer  into  the  equation  of 
state,  writing  the  resultant  expression  in  differential  form,  and  then 
finally  equating  that  expression  to  the  hydrostatic  balance  equation,  Eq 
(A-4) ,  to  obtain  the  following  ordinary  differential  equation  for  the 
density, 

dp  «  -g  dhp  where  C  ■  RO2  TSL  (A-7) 

P  C 

This  equation  was  integrated  between  hp  =  11  km  and  20  km  to  determine 
the  variation  of  Layer-2  density  as  a  function  of  geopotential  altitude 
which  is  given  as  follows  in  terms  of  the  difference  in  altitude  above 
the  tropopause,  (hp  -  h  ), 

p  ■  p*  exp  [  -(g/C)  (hp  -  h*)]  (A-8) 
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Non-dimensionalization  of  the  above  expression  is  accomplished  by 
dividing  through  by  the  sea  level  density  (  ), 

a  £  -  P  «  P  *  exp  [  -(g/C)  (hp  -  h*)]  (A-9) 

^SL  ^SL 

A 

However,  a  numerical  value  of  P  must  still  be  obtained  to  define  the 
density  variation  with  geopotential  altitude.  Fortunately,  at  this 
point  both  the  non-dimensional  pressure  and  temperature  variation  with 
geopotential  altitude  are  known  in  terms  of  basic  parameters. 
Additionally,  since  the  equation  of  state  is  expressible  in  terms  of  the 
non-dimensional  quantities, 

6  *  o  •  0  (A-10) 


the  value  of  P 


* 


can  be  explicitly  obtained  in  the  following  manner: 


“  P 


SL 


.  0  4.2559121 


(A-li) 


Consequently,  Eq  (A— 9)  now  can  be  expressed  in  its  final  form  as 
follows: 


o2  «  0.29707289  exp  [ -4. 8063758(10)"5  (hp  -  h*)]  (A-12) 
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+A/m*t*rj?n*  ^»jv>vyg.r>s<  'Uj-«<.^  ^  m-.^-^,,-.^,-^-^  - 


vrqF*^Gm*e'Vj\*4tct;$& *psr^  '?V$ 


1  : 


Since  the  non-dimensional  temperature  and  density  variation  with 
geopotential  altitude  for  the  second  layer  are  now  known,  the 
corresponding  expression  for  the  non-dimensional  pressure  variation  with 
altitude  can  be  determined  from  the  non-dimensional  equation  of  etate, 
Eq  (A-10) ,  as  indicated  by. 


62  “  02  02 


0.22335881  exp  [-4.8063758(10)~5  (h  -  h*)] 


(A-13) 


Similarly,  the  non-dimensional  density  variation  of  the  first  layer 
is  obtainable  from  the  non-dimensional  equation  of  state  as  follows: 


61  -  q^.2559121 


(A-14) 


A.  3  Summary 

The  following  expressions  summarize  the  non-dimensional  variation 
of  temperature,  pressure,  and  density  within  the  first  two  layers  of  the 
1962  United  States  standard  atmosphere. 


Layer- 1 


0-  1  -  6.8755856(10)“6  h_ 


5.2559121 


4.2559121 
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0 


0.75186535 


as 

5  -  0.22335881  exp  [ -4.8063758(10)“5  (hp  -  h* 
°  “  0.29707289  exp  [ -4.8063/58(10)~5  (hp  -  h* 


1S3 


Appendix  1$ 

Gaussian  Jet  Velocity  Profile 

In  addition  to  the  uniform  ("top-hat")  velocity  profile 
accommodated  within  the  computer  simulation,  a  truncated  bivariate 
normal  distribution  was  also  included  for  possible  use  to  approximate 
the  exhaust  gas  velocity  field  as  it  was  forcibly  being  ejected  through 

a  rectangular  diffuser  exit  port  into  a  high  subsonic  crossflow.  The 

truncated  bivariate  density  function  was  expressed  in  the  following  form 
(Ref  B.l) : 

,  „  2  2 
K  exp  f  -1  ,x-^  '  + 

2  to  a  1  2  (  o  ) 

x  y  v.  L'  x 

for  |x|  &  a,  |y|  <  b;  (B-l) 

0  ,  for  |x|  >  a,  jy|  >  b  . 

where  :  x^ ,  are  the  mean  values  of  x  and  y 

respectively, 

<J  ,  o  are  the  standard  deviations  of  x  and  y 
x  y 

respectively,  and, 

2a,  2b  are  the  dimensions  of  the  rectangular  jet 
orifice,  2a  being  aligned  with  the  freestream 
direction. 
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The  parameter  K  in  Eq  (B-l)  accounts  for  the  fact  that  the  orifice  is  of 
finite  dimension  and  that  the  approximation  of  the  jet  injection 
velocity  by  a  Gaussian  distribution  over  this  finite  region  requires  a 
truncation  factor  to  compensate  for  the  neglected  part  of  the  continuous 
Gaussian  density  function.  The  value  of  K  is  calculated  from  the 
following  constraint  equation, 


f(x,y)  dx  dy 


1 


(B-2) 


When  the  expression  of  the  density  function,  Eq  (B-l),  is  substituted 
into  the  above  constraint  equation,  and,  the  indicated  integrations  are 
carried  out,  the  following  equation  results: 

a  b 

f(x,y)  dx  dy  = 


/k  erff/ip  x+yTa]-  K  erf[~/rpx  -  x/Ta]^  erf  [v^p  y  +V/Fb"|  + 

\  L  2  0  .  L  2  a  .  /  L  2  a 

'  x  x  /  y 

K  erf  pTp  x  -  v/F  al  -  K  erf  1^2  y  x  +  /2_a'l\  erff/Iy  y  -/Tb]  (B-3) 

2ax  J  L  2  a x  J/  L  2oy 

Given  a  rectangular  jet  of  dimension  (2a)  x  (2b),  if  the  mean  values, 
and  y^,  are  assumed  to  be  identically  zero,  then  the  following 
expression  results  from  which  the  required  value  of  K  can  be  obtained 
directly,  namely. 
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K  erf J  a  1  erf  b  -  1  (B-4) 

T*5! !  Wy 

It  is  easily  verified  that  if  the  limits  in  Eq  (B-3)  are  allowed  to 
approach  those  associated  with  the  true  bivariate  Gaussian  distribution, 
that  is,  if  both  a  and  be  were  allowed  to  become  infinitely  large  in  the 
limit,  and,  if  both  standard  deviations  were  constrained  to  remain 
positive,  then  the  integral  simplifies  as  required  to  the  following 
identity,  that  is, 

K  ■=  1  (B-5) 

If  each  half  dimension  of  the  rectangular  orifice  is  assumed  to 
represent  the  standard  deviation  in  the  respective  coordinate  direction, 
that  is,  if 

a  =  0  and  b  «=  o  (B-6) 

x  y 

and  furthermore,  if  the  mean  values  of  x  and  y  are  coincident  with  the 
values  associated  with  the  geometric  center  of  the  orifice,  then  the 
truncation  correction  factor  for  the  truncated  bivariate  density 
function  is  expressible  as 

K  -  1  (B-7) 

erf (a)  erf(S) 

where  a  -yfljl  •  S,  from  which  the  value  of  K  for  this  special  case 
is  immediately  determined, 
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■fj 

ft-' 

I 

& 


£ 


K 


2.1456238 


(B-8) 


The  truncated  bivariate  Gaussian  velocity  distribution  associated  with 
the  assumptions  expressed  in  Eq  (B-6)  possesses  a  relatively  flat 
profile  which  is  characteristic  of  fully  developed  turbulent  jet3. 
Consequently,  Eq  (B-l)  with  the  value  of  K  expressed  in  Eq  (B-8)  was 
utilized  to  define  the  "normal"  form  of  the  velocity  profile  of  the  jet 
within  the  computer  simulation. 


Appendix  0 

Gas  Transport  Properties 


Modeling  real  gases  at  elevated  temperatures  requires  the 
additional  complexity  incurred  by  determining  the  gas  mixture  transport 
properties  of  viscosity  (P  ),  thermal  conductivity  (A.  ),  and  molecular 
diffusion  coefficient  (P).  The  corresponding  properties  of  the  pure 
component  gases  were  approximated  by  using  the  well  established  Lennard- 
Jones  12-6  intermolecular  potential  function  (Ref  C.l)  which  is  given  by 


<j)(r)  *  4e 

a 

12 

a 

6 

r 

r 

which  is  clearly  a  function  of  only  two  parameters  for  any  gi/en  gas;£  , 
the  maximum  energy  of  attraction  (well  depth  of  the  potential  function) 
expressed  in  ergs,  and  O,  the  zero  energy  collision  diameter  expressed 
in  angstroms.  Svehla  (Ref  C.2)  documented  these  force  constants  (£  /k 
and  a  where  k  is  the  Boltzmann  constant)  for  a  large  number  of  gaseous 
species.  Additional  values  can  be  obtained  from  the  Hirschf elder , 
Curtiss  and  Bird  text  (Ref  C.l)  for  the  more  ordinary  gases  and  from  the 
report  by  Liley  (Ref  C.3)  for  the  more  exotic  species  more  directly 
associated  with  the  exhausts  of  chemical  lasers.  The  Lennard-Jones 
force  constants  utilized  within  the  computer  simulation  are  tabulated  in 

_  (1»D* 

Table  C.l.  Tables  of  the  required  collision  integrals,  S2  and 

(2,2)* 

ft  ,  used  in  the  requisite  expressions  for  gas  transport  properties 
were  obtained  from  Ref  (C.l)  despite  some  errors  allegedly  inherent  in 
them  as  reported  by  Liley  (Ref  C.3).  Empirical  expressions  for  these 
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integrals  whic'..  are  well  suited  to  numerical  programming  are  presented 
in  Refs  C.4  and  C.5;  however,  these  expressions  neither  are  as  accurate 
as  the  tabular  values  nor  are  they  applicable  over  the  entire  range  of 
values  encountered  during  the  numerical  simulation.  Consequently,  the 
tabular  values  of  Hirschf elder,  Curtiss  and  Bird  were  retained. 


C.l  Mixture  Viscosity 

The  first  approximation  of  the  molecular  viscosity  for  a  pure  gas 
may  be  written  (Ref  C.l,  Eq  8.2-18)  as, 


M  -  266-93  (iO)-7 


c^2-2** 


,  poises 


(C— 2 ) 


A2,2)* 

where  M  is  the  molecular  weight  of  the  gas  species,  and  M  is  a 
tabulated  function  of  the  reduced  temperature,  T*  =  T/( e  /k  )  ,  where 
T  is  absolute  temperature  expressed  in  degrees  Kelvin. 

Following  the  procedure  outlined  in  Ref  C.l,  the  updated  value  of 
the  pure  gas  component  viscosity  was  obtained  by  using  Eq  8.2-19  of  that 
reference  which  is  given  by, 


y*  V 


(C-  ’) 


where  the  function  f^K'  is  the  kc  -order  correction  factor  whose  value 
is  approximately  equal  to  one  (1.).  Values  of  this  correction  are 
tabulated  (Table  I-P)  as  functions  of  T*  in  the  text  of  Hirschf  elder , 


Curtiss  and  Bird. 


The  rigorously  derived  expressions  according  to  the  Chapman-Ens kog 


theory  for  the  viscosity  of  a  multi- component  gas  mixture  are  presented 
in  Ref  C.l,  Eqs  8.2-25  -  28;  however,  they  involve  the  ratios  of 
determinants  whose  elements  are  neither  easily  nor  efficiently 
programmable  for  gas  mixtures  comprised  of  more  than  two  distinct 
species.  The  rigorous  expansion  can  be  closely  approximated  for  non¬ 
polar  gases  at  low  pressures  by  the  following  series  (Ref  C.5,  Eq  9- 
5.1), 


y  mix  =  X  ^i  (C"4) 

i=l  N 

Z  xj  *ij 
j-1 


where  are  the  component  pure  gas  viscosities,  and,  x^  are  the  mole 
fractions  of  each  species. 

Furthermore,  utilization  of  Wilke's  approximation  for  the  parameter 
(Ref  C.5,  Eq  9-35)  as  given  by. 


where  are  the  component  molecular  weights,  results  in  an  easily 
programmable  expression  for  the  viscosity  of  the  gaseous  mixture. 
Brokaw  (Ref  C.6)  developed  alignment  charts  for  both  the  component 
viscosities  and  the  expression  given  by  Eq  (C-5).  A  method  for 
calculating  the  viscosity  of  a  gas  mixture,  some  components  of  which  are 
polar,  was  detailed  by  Brokaw  in  Ref  C.7. 
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C.2  Mixture  Thermal  Conductivity 

The  first  approximation  for  the  thermal  conductivity  of  a  pure 
monatomic  gas  may  be  written  in  the  following  form  (Ref  C.l,  Eq  8.2-31): 


X  «  1989.1  (10)"'  _/HT  #  cal/(cm-sec-°K) 


(C-6) 


a2  ft 


(2,2)* 


-(2,2)* 

again  where  M  is  the  pure  gas  species  molecular  weight,  ft  is  the 
collision  integral  tabulated  as  a  function  of  the  reduced  temperature, 
T  ,  and  T  is  expressed  in  degrees  Kelvin.  Updating  this  first 
approximation  to  the  pure  gas  thermal  conductivity,  the  methodology  of 
Hirschf elder,  Curtiss,  and  Bird  was  utilized  (Ref  C*l,  Eq  8. 2 '32)  to 
obtain  the  order  approximation  given  by, 


[X  ]k  -  X  *  f(Ak)  (C-7) 

where  f ^  is  the  ktk  order  correction  factor  tabulated  in  Ref  C.l, Table 
I-P,  as  a  function  of  T*  and  it  too  is  approximately  equal  to  unity. 

Analagously  to  the  mixture  viscosity  calculations,  the  rigorous 
application  of  Chapman-Enskog  theory  for  non-polar  gases  results  in  the 
mixture  thermal  conductivity  being  expressed  again  as  a  ratio  of 
determinants  (Ref  C.l,  Eq  8.2-43)  whose  elements  are  equally  unwieldy 
for  computational  purposes.  Instead,  the  Cheung,  Bromley,  and  Wilke 
method  (Ref  C.5,  Eq  10-6.1)  utilizing  the  familiar  Wassiljewa  form  of 
the  approximation  was  employed  to  obtain  the  thermal  conductivity  for 
the  gas  mixture  composed  of  non-polar  gases,  that  is, 
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^L-*r*^,^<n>-,St;^,V'4fejn  -*■ 


X  mix  "  2^  XiX  i 
i»l  N 


(C-8) 


2  Xj  Aij 
j“l 

where  ^  ^  is  the  thermal  conductivity  of  the  ic^  species  given  by  Eq  (C- 
6)  ,  x^  and  x^  are  the  mole  fractions  of  species  "i"  and  "j" 
respectively,  and  A^j  is  a  parameter  which  is  temporarily  unspecified. 

The  method  proceeded  by  splitting  the  thermal  conductivity  of  each 
species  into  two  distinct  parts. 


X  _  X  *  .  X  ** 

i  i  +  i 


(C-9) 


where  the  first  part,  X  ,  accounts  for  the  contribution  of  the 
monatomic,  or  translational,  thermal  conductivity  of  the  species,  and 
the  second  part,  X  ^  ,  accouunts  for  the  contribution  of  the  polyatomic 
thermal  conductivity  by  internal  energy  diffusional  transport , or 
internal  degrees  of  freedom  (Ref  C.2,  Appendix  B).  After  applying  the 
Euken  corrections  to  ^  to  account  for  the  presence  of  polyatomic 
gases,  the  following  expressions  lor  this  quantity  resulted  (Ref  C.4,  Eq 
10-40) : 


f  K  ^  ,  monatomic 

X  i  “  '  X  i  j1/  [  1  +  0.35  (a?  /  R  -  2^j  |  ,  linear  (C-10) 

Xi  |l/[l  +0-25  (c  /  R  -  2)]  |  ,  non-linear 


"k'k 

After  X  ^  was  obtained  from  Eq  ( C-9)  and  the  unspecified  parameter  A^ 
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Eq  (C-5),  the  mixture 


was  assumed  to  be  well  represented  by  <j>  ^ ^  , 

thermal  conductivity  was  expanded  in  the  following  manner: 


X  mix  =  £X  1  xi 


i=l  N 


0.125 


E  <vmi> 


+  Xx  r  xi  (c_ii) 

i=l  ii 

E  «, 

j=l 


where 


My  =  (  M±  +  Mj  )  /  2 


( C— 12) 


C . 3  Mixture  Molecular  Diffusion  Coefficient 

The  first  approximation  of  the  coefficient  of  diffusion  for  a 
binary  gas  mixture  comprised  of  species  "i"  and  "j"  is  given  by  the 
following  expression  (Ref  C.l,  Eq  8.2-44), 


-  0.0026280  \/T3  (M1-mi)/(2MiM1)  ,  (cm2/sec) 

(1,D* 


P  aij 


(C-13) 


where  P  is  the  pressure  in  atmospheres,  T  is  the  temperature  in  degrees 

(1,1)* 

Kelvin,  is  the  collision  integral  evaluated  at  the  temperature 


Lij] 


(C-14) 


where 


is  the  pair  molecular  potential  energy  parameter  given  by, 
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(C-15) 


°ij  -  (ci-i-aj)/2  ,  (angstroms) 

this  approximation  of  V  ^  was  similarly  updated  by  multiplication  with 
the  k*1*1  order  correction  factor*  fffit  given  ay  (Ref  C.l,  Eq  8.2-45) 

t»i j  4  ■  '  f!pk)  <c-16> 

where  f^  is  also  tabulated  in  Table  I-P  of  Ref  C.l  as  a  function  of 
the  :.ad»ced  temperature,  T  .  This  correction  factor  also  is 
approxL-atedly  equal  to  unity. 

Calculation  of  the  diffusion  coefficient  for  a  mixture  of  gases  is 
extremely  difficult  if  not  nearly  impossible  for  gaseous  mixtures 
representative  of  exhausts  from  chemical  laser  systems.  However,  if  one 
starts  with  the  Stef an-Maxwell  equations  given  by  (Ref  C.8,  Eq  18.4-19), 


SB 


N 


Vi  } 


(C-l? ) 


where  x.^  is  the  mole  fraction  o£  species  "i",  is  the  binary  mixture 
diffusion  coefficient  for  the  i-j  species  gaseous  pair,  and  v,  and  Sj 
are  the  velocities  of  species  "j"  and  "i"  respectively,  and  if,  it  ts 
further  assumed  that  one  homogeneous  gas  diffuses  into  another 
homogeneous  gas  (for  example,  air  diffusing  into  the  exhaust  jet  gases), 
then  it  becomes  possible  to  derive  the  diffusion  coefficient  of  the 
binary  mixture  comprised  of  both  the  homogeneous  gases.  The  resulting 
expression  of  the  molecular  diffusion  coefficient  for  this  binary 
gaseous  pair  is  then  given  by  (Ref  C.8,  Eq  18.4-25), 
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(C-18) 


(cm^/sec) 


where  x ^  and  x^  are  the  mole  fractions  of  the  diffusing  gas  and  the 
exhaust  gases  respectively,  and  V  ^  is  the  binary  mixture  diffusivity 
of  the  diffusing  gas  with  each  of  the  exhaust  species  gases  comprising 
the  laser  exhaust  which  is  given  for  each  cf  these  pairs  by  Eq  (C-13). 


C . 4  Unit  Conversion 

The  units  utilized  within  the  computer  simulation  were 
predominantly  English  engineering  units.  The  units  associated  with  the 
transport  properties  of  the  gases  were  initially  computed  in  various 
forms  of  the  Metric  system  thereby  necessitating  numerical  conversion. 
The  conversion  factor  compendium  of  Mechtley  (Ref  C.9)  was  used  to 
facilitate  conversion  of  the  derived  properties  to  compatible  English 
engineering  unite.  The  conversion  tables  in  Ref  C.9  are  equally  usable, 
however,  the  corresponding  factors  are  not  carried  to  the  same  precision 
es  those  provided  by  Mechtly. 


Table  C.l 

Lennard  Jones  (12-6)  Force  Constants 


>ecies 

Molecular 

e/k 

0 

0 

Weight 

(°K) 

(A) 

CF4 

83.00475 

134.0 

4.662 

He 

4.00260 

10.22 

2.551 

°2 

4.0028204 

35.20 

2.952 

HF 

20.006303 

330.0 

3.418 

DF 

21.0125 

199.1 

2.826 

n2 

28.016 

71.4 

3.798 

Air 

28.9644 

78.6 

3.711 
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Appendix  D 

Thermodyramic  Properties 


Gas  mixture  thermodynamic  properties  were  calculated  from  component 
species  thermodynamic  properties  based  upon  the  respective  mass 
fractions  of  each  species.  Values  of  specific  heat  at  constant 
pressure,  Cp,  and  sensible  enthalpy,  ^“^29 8. 15*  were  tabulated  for  each 
gas  mixture  component  as  functions  of  temperature  from  200  -  2500°K  . 
The  values  for  the  exhaust  gas  species  were  obtained  from  the  JANAF 
Thermochemical  Tables  (Ref  D.l).  The  corresponding  values  for  air  were 
obtained  from  two  separate  sources;  sensible  enthalpies  were  obtained 
from  the  tables  of  Keenan  and  Kaye  (Ref  D.2)  while  specific  heats  were 
obtained  from  the  text  by  Eckert  and  Drake  (Ref  D.3,  Table  B-4).  Tables 
D.l  and  D.2  list  the  values  of  cp  and  for  every  gas  component 
required  by  the  computer  simulation.  Required  unit  conversions  were 
again  accomplished  using  the  compendium  of  conversion  factors  compiled 


by  Mechtly  (Ref  C.9). 
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